forked from bpartopour/PBG
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathradial_voidage.py
170 lines (154 loc) · 6.3 KB
/
radial_voidage.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
import bpy
from mathutils import Vector
from mathutils.bvhtree import BVHTree
import numpy as np
import math
import os
import parameters
###### function q_inside checks if a point is inside the "solid" packing based on ray-casting algorithm ###
###### ray casting algorithm can be mathematically proved by Jordan curve theorm #########
def q_inside(obj, point_in_object_space):
direction = Vector((1,0,0))
epsilon = direction * 1e-6
count = 0
result, point_in_object_space, normal, index = obj.ray_cast(point_in_object_space, direction)
while result:
count += 1
result, point_in_object_space, normal, index = obj.ray_cast(point_in_object_space + epsilon, direction)
return (count % 2) == 1
def radial_voidage():
bpy.ops.object.delete()
bpy.ops.import_mesh.stl(filepath= parameters.file_path, axis_forward='Y', axis_up='Z', filter_glob="*.stl", global_scale=1.0, use_scene_unit=True, use_facet_normal=False)
for ob in bpy.context.scene.objects:
if ob.type == 'MESH':
ob.select = True
bpy.context.scene.objects.active = ob
for obj in bpy.data.objects:
obj.name = 'bed'
# renaming the stl file
obj = bpy.data.objects["bed"]
# changing to edit mode to make all the face normals consistant and facing outwards
bpy.ops.object.editmode_toggle()
bpy.ops.mesh.normals_make_consistent(inside = False)
bpy.ops.object.editmode_toggle()
#setting the origin to the geometry and to the space origin
bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')
bpy.context.object.location = 0,0,0
p_l = obj.matrix_world.translation
u = bpy.context.object.dimensions
print("The bed dimensions (X,Y,Z) are:", u[0], u[1], u[2])
R_tube = (max(u[0],u[1]))/2.0
d = u[2] - 8*parameters.particle_radius
if d <= (2*parameters.particle_radius):
decision = ''
while decision != 'Yes' and decision != 'No':
decision = input('WARNING: The bed is not long enough for accurate radial voidage calculation, do you want to continue? Yes, No?')
if decision != 'Yes' and decision != 'No':
print("your answer is not acceptable, please enter Yes or No")
if decision == 'Yes':
d = u[2]
r_p = [0]*99
w_d = [0]*99
j=1
dp= 2*parameters.particle_radius
# setting max and min height of the radial surface
Lmin = -d/2
Lmax = d/2
# generating a cylindrical surface mesh
w = 2*math.pi
r_p = [0]*99
w_d = [0]*99
for i in range(0,99):
#print('loop number is', i)
R = j*R_tube
step = 2*math.asin(dp/20/2/R)
list_t = [t for t in np.arange(0,w,step)]
L_face = dp/20
#print(L_face)
Area_face = pow(L_face,2)
list_x_y = [(R*math.cos(x_y),R*math.sin(x_y)) for x_y in list_t]
list_z = [z for z in np.arange(Lmin,Lmax,L_face)]
point_co =[(x,y,z) for (x,y) in list_x_y for z in list_z]
point_nu = len(point_co)
count = 0
#print('number of points',point_nu)
d = False
for p in range(0,point_nu):
point_in_object_space = Vector((point_co[p]))
d = q_inside(obj, point_in_object_space)
if d:
count += 1
ring_area = Area_face * point_nu
packed_area = Area_face * count
r_p[i] = 1-(count/point_nu)
#print('eps', r_p[i])
w_d[i] = j
j = j - 0.01
if i%10 == 0:
print((i),'% progress...')
print('100% progress...Radial voidage calculation is finished!')
print('writing the data...')
w_d[0] = 1.0
Distance = [(R_tube*(1-x))/(2*parameters.particle_radius) for x in w_d]
file_name = filepath=parameters.file_path_2
f=open(file_name,'w')
for k in range(99):
d = Distance[k]
z = r_p[k]
f.write("Distance %.2f Voidage %.2f \n"%(d, z))
f.close()
# print('Done! :)')
else:
print('Radial voidage calculation has stopped, goodbye!')
elif d > (2*parameters.particle_radius):
r_p = [0]*99
w_d = [0]*99
j=1
dp= 2*parameters.particle_radius
# setting max and min height of the radial surface
Lmin = -d/2
Lmax = d/2
# generating a cylindrical surface mesh
w = 2*math.pi
r_p = [0]*99
w_d = [0]*99
for i in range(0,99):
#print('loop number is', i)
R = j*R_tube
step = 2*math.asin(dp/20/2/R)
list_t = [t for t in np.arange(0,w,step)]
L_face = dp/20
#print(L_face)
Area_face = pow(L_face,2)
list_x_y = [(R*math.cos(x_y),R*math.sin(x_y)) for x_y in list_t]
list_z = [z for z in np.arange(Lmin,Lmax,L_face)]
point_co =[(x,y,z) for (x,y) in list_x_y for z in list_z]
point_nu = len(point_co)
count = 0
#print('number of points',point_nu)
d = False
for p in range(0,point_nu):
point_in_object_space = Vector((point_co[p]))
d = q_inside(obj, point_in_object_space)
if d:
count += 1
ring_area = Area_face * point_nu
packed_area = Area_face * count
r_p[i] = 1-(count/point_nu)
#print('eps', r_p[i])
w_d[i] = j
j = j - 0.01
if i%10 == 0:
print((i),'% progress...')
print('100% progress...Radial voidage calculation is finished!')
print('writing the data...')
w_d[0] = 1.0
Distance = [(R_tube*(1-x))/(2*parameters.particle_radius) for x in w_d]
file_name = filepath=parameters.file_path_2
f=open(file_name,'w')
for k in range(99):
d = Distance[k]
z = r_p[k]
f.write("Distance %.2f Voidage %.2f \n"%(d, z))
f.close()
# print('Done! :)')