Skip to content

Commit fe543d5

Browse files
chyalexchengchyalexcheng
committed
update GetNetWorkData.py
1 parent 1aa6443 commit fe543d5

File tree

1 file changed

+117
-84
lines changed

1 file changed

+117
-84
lines changed

‎data_processing/GetNetworkData.py

Lines changed: 117 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -10,73 +10,81 @@
1010
import os, glob
1111
import h5py
1212

13-
contact_keys = [
14-
'O.materials[0].young', # young's modulus = 10^E
15-
'O.materials[0].poisson', # poisson's ratio
16-
'O.materials[0].frictionAngle', # sliding friction
13+
contact_keys = {
14+
'young_modulus':'O.materials[0].young', # young's modulus = 10^E
15+
'possion_ratio':'O.materials[0].poisson', # poisson's ratio
16+
'friction_angle':'O.materials[0].frictionAngle', # sliding friction
17+
'particle_density':'O.materials[0].density', # density
18+
}
19+
20+
particle_keys = {
21+
'radius':'shape.radius', # particle radius
22+
# 'mass':'shape.mass', # particle mass can be computed as 4/3*density*pi*radius**3
23+
}
24+
25+
sampling_keys = [
26+
'pressure', # confining pressure
27+
'initial_friction', # initial friction (for generating microstructure with different void ratio)
28+
'compressive_strain_rate', # rate of compressive strain d{\epsilon_p}
29+
'shear_strain_rate', # rate of shear strain d{\epsilon_q}
1730
]
1831

1932
## The following list contains macro-variables. For drained triaxial, the macro-input and -output are like below. Note that for the initial load step, all macro-variables are input.
20-
macro_keys = [
33+
macro_keys = {
2134
## macroscopic input (strain can be computed from log(l_x^i/l_x^0) and void ratio from the sum of particle volume)
22-
'O.cell.hSize[2,2]', # domain size in z
23-
'getStress()[0,0]', # stress in x
24-
'getStress()[1,1]', # stress in y
35+
'domain_size_z':'O.cell.hSize[2,2]', # domain size in z
36+
'stress_x':'getStress()[0,0]', # stress in x
37+
'stress_y':'getStress()[1,1]', # stress in y
2538
## macroscopic output
26-
'O.cell.hSize[0,0]', # domain size in x
27-
'O.cell.hSize[1,1]', # domain size in y
28-
'getStress()[2,2]', # stress in z
39+
'domain_size_x':'O.cell.hSize[0,0]', # domain size in x
40+
'domain_size_y':'O.cell.hSize[1,1]', # domain size in y
41+
'stress_z':'getStress()[2,2]', # stress in z
2942
## others
30-
'O.dt', # time increment
31-
'O.iter', # number of iterations
32-
]
43+
'time':'O.time', # time increment
44+
'num_time_steps':'O.iter', # number of time steps
45+
}
3346

34-
micro_keys_bodies = [
47+
micro_keys_bodies = {
3548
## particle info (in b.shape and b.state)
36-
'id',
37-
'shape.radius',
38-
'state.mass',
39-
'state.pos[0]',
40-
'state.pos[1]',
41-
'state.pos[2]',
42-
'state.vel[0]',
43-
'state.vel[1]',
44-
'state.vel[2]',
45-
'state.angVel[0]',
46-
'state.angVel[1]',
47-
'state.angVel[2]',
48-
'state.refPos[0]', # position at t=0, with respect to the origin
49-
'state.refPos[1]',
50-
'state.refPos[2]',
51-
]
52-
micro_keys_inters = [
49+
'position_x':'state.pos[0]',
50+
'position_y':'state.pos[1]',
51+
'position_z':'state.pos[2]',
52+
'velocity_x':'state.vel[0]',
53+
'velocity_y':'state.vel[1]',
54+
'velocity_z':'state.vel[2]',
55+
'angular_velocity_x':'state.angVel[0]',
56+
'angular_velocity_y':'state.angVel[1]',
57+
'angular_velocity_z':'state.angVel[2]',
58+
}
59+
60+
micro_keys_inters = {
5361
## interaction connectivity info (in inter)
54-
'id1',
55-
'id2',
62+
'contact_pair_ID1':'id1',
63+
'contact_pair_ID2':'id2',
5664
## interaction geometry info (in inter.geom)
57-
'phys.ks', # tangential stiffness
58-
'phys.kn', # normal stiffness
59-
'geom.penetrationDepth', # overlap between spheres
60-
'geom.shearInc[0]', # shear increment x between particles
61-
'geom.shearInc[1]', # shear increment y between particles
62-
'geom.shearInc[2]', # shear increment z between particles
63-
'geom.contactPoint[0]', # x, y, z, in the cross section of the overlap
64-
'geom.contactPoint[1]',
65-
'geom.contactPoint[2]',
65+
'tangential_stiffness':'phys.ks', # tangential stiffness
66+
'tangential_stiffness':'phys.kn', # normal stiffness
67+
'overlap':'geom.penetrationDepth', # overlap between spheres
68+
'shear_increment_x':'geom.shearInc[0]', # shear increment x between particles
69+
'shear_increment_y':'geom.shearInc[1]', # shear increment y between particles
70+
'shear_increment_z':'geom.shearInc[2]', # shear increment z between particles
71+
'contact_point_x':'geom.contactPoint[0]', # x, y, z, in the cross section of the overlap
72+
'contact_point_y':'geom.contactPoint[1]',
73+
'contact_point_z':'geom.contactPoint[2]',
6674
## interaction physics info (in inter.phys)
67-
'phys.usElastic[0]', # elastic component of the shear displacement x
68-
'phys.usElastic[1]', # elastic component of the shear displacement y
69-
'phys.usElastic[2]', # elastic component of the shear displacement z
70-
'phys.usTotal[0]', # total shear displacement x
71-
'phys.usTotal[1]', # total shear displacement y
72-
'phys.usTotal[2]', # total shear displacement z
73-
'phys.shearForce[0]', # shear foce x
74-
'phys.shearForce[1]', # shear foce y
75-
'phys.shearForce[2]', # shear foce z
76-
'phys.normalForce[0]', # normal force x
77-
'phys.normalForce[1]', # normal force y
78-
'phys.normalForce[2]', # normal force z
79-
]
75+
'elastic_shear_x':'phys.usElastic[0]', # elastic component of the shear displacement x
76+
'elastic_shear_y':'phys.usElastic[1]', # elastic component of the shear displacement y
77+
'elastic_shear_z':'phys.usElastic[2]', # elastic component of the shear displacement z
78+
'total_shear_x':'phys.usTotal[0]', # total shear displacement x
79+
'total_shear_y':'phys.usTotal[1]', # total shear displacement y
80+
'total_shear_z':'phys.usTotal[2]', # total shear displacement z
81+
'shear_force_x':'phys.shearForce[0]', # shear foce x
82+
'shear_force_y':'phys.shearForce[1]', # shear foce y
83+
'shear_force_z':'phys.shearForce[2]', # shear foce z
84+
'normal_force_x':'phys.normalForce[0]', # normal force x
85+
'normal_force_y':'phys.normalForce[1]', # normal force y
86+
'normal_force_z':'phys.normalForce[2]', # normal force z
87+
}
8088

8189
unused_keys_sequence = [
8290
]
@@ -89,7 +97,10 @@
8997
DATA_DIR = '/home/cheng/DataGen/'
9098
STORE_EDGE_FEATURES = False
9199

92-
100+
# load the sampling variables
101+
state_sampling = [log10(0.1e6), log10(0.01)]
102+
path_sampling = [0, 0]
103+
93104
def main(pressure, experiment_type,numParticles):
94105
data_dir = DATA_DIR + f'{pressure}/{experiment_type}/{numParticles}/'
95106
if not os.listdir(data_dir):
@@ -99,62 +110,83 @@ def main(pressure, experiment_type,numParticles):
99110
# get DEM state file names
100111
simStateName = 'simState_' + experiment_type + f'_*_{numParticles}_0.yade.gz'
101112
file_names = glob.glob(data_dir + '/' + f'{simStateName}')
113+
102114
# get the list of sample IDs
103115
samples = [int(f.split('.yade.gz')[0].split('_')[-3]) for f in file_names]
104116
simStateName = data_dir + '/' + 'simState_' + experiment_type
117+
print(f'Number of samples is {len(samples)}')
105118

106119
# name the HDF5 file
107-
h5file_name = data_dir + '/' + 'simState_' + experiment_type + f'_all_{numParticles}' + '_graphs.hdf5'
120+
h5file_name = data_dir + '/' + 'simState_' + experiment_type + f'_all_{numParticles}' + '_graphsClean.hdf5'
108121
if os.path.exists(h5file_name):
109122
os.remove(h5file_name)
110123
f_all = h5py.File(h5file_name, 'a')
111124

112-
# store the keys and constant parameter and
113-
f_all['contact_keys'] = contact_keys
114-
f_all['macro_keys'] = macro_keys
115-
f_all['micro_keys_bodies'] = micro_keys_bodies
116-
f_all['micro_keys_inters'] = micro_keys_inters
125+
# store the variable keys
126+
f_all.attrs['contact_keys'] = list(contact_keys.keys())
127+
f_all.attrs['radius'] = list(particle_keys.keys())
128+
f_all.attrs['sampling_variables'] = sampling_keys
129+
f_all.attrs['macro_input_features'] = [list(macro_keys.keys())[i] for i in [0,1,2]]
130+
f_all.attrs['macro_output_features'] = [list(macro_keys.keys())[i] for i in [3,4,5]]
131+
f_all.attrs['time'] = list(macro_keys.keys())[6:]
132+
f_all.attrs['node_features'] = list(micro_keys_bodies.keys())
133+
f_all.attrs['edge_features'] = list(micro_keys_inters.keys())
134+
135+
# load a simulation state to get the contact parameters
117136
O.load(f'{simStateName}_1_{numParticles}_1.yade.gz')
118-
contact_tensor = np.array([float(eval(key)) for key in contact_keys])
119-
f_all['contact_params'] = contact_tensor
137+
contact_tensor = np.array([float(eval(key)) for key in contact_keys.values()])
138+
139+
# store metadata universal to all samples
140+
f_all_meta = f_all.create_group('metadata')
141+
f_all_meta['contact_params'] = contact_tensor
142+
f_all_meta['num_nodes'] = int(numParticles)
143+
f_all_meta['mean_radius'] = 0.5
144+
f_all_meta['dispersion_radius'] = 0.4
120145

121146
# load YADE and store data in f_all
122-
for sample in samples:
123-
simStateName = data_dir + '/' + 'simState_' + experiment_type
147+
for sample in samples[:2]:
124148
# load the initial time step
125149
try:
126150
O.load(f'{simStateName}_{sample}_{numParticles}_0.yade.gz'); O.step()
127-
except IOError:
128-
print('IOError', f, pressure)
151+
except RuntimeError:
152+
print(f'RuntimeError encountered for {simStateName}_{sample}_{numParticles}_0.yade.gz')
129153
continue
130154

131155
## get DEM state file names
132156
file_names = glob.glob(f'{simStateName}_{sample}_{numParticles}_*.yade.gz')
133157
# get the list of loadstep IDs
134158
steps = sorted([int(f.split('.yade.gz')[0].split('_')[-1]) for f in file_names])
135-
159+
print(f'Steps: {steps}')
160+
161+
# create a group per sample
136162
f_sample = f_all.create_group(f'{sample}')
137-
f_sample['num_steps'] = len(steps)
138-
139-
for step in steps[1:]:
163+
# store metadata specific to each sample
164+
f_sample_meta = f_sample.create_group('metadata')
165+
for m,key in enumerate(sampling_keys): f_sample_meta[key] = (state_sampling+path_sampling)[m]
166+
f_sample_meta['num_steps'] = len(steps)
167+
f_sample_meta['radius'] = np.array([float(eval('b.'+const_body_key)) for b in O.bodies for const_body_key in particle_keys.values()]).reshape([int(numParticles),1])
168+
# store the time sequence of each sample
169+
f_sample_time = f_sample.create_group('time_sequence')
170+
# loop over time
171+
for i,step in enumerate(steps[:2]):
140172
# load YADE at a given time step
141173
try:
142174
O.load(f'{simStateName}_{sample}_{numParticles}_{step}.yade.gz'); O.step()
143-
except IOError:
144-
print('IOError', f, pressure)
175+
except RuntimeError:
176+
print(f'RuntimeError encountered for {simStateName}_{sample}_{numParticles}_{step}.yade.gz')
145177
continue
146178
## macroscopic data (domain size, stress, etc)
147-
macro_tensor = np.array([float(eval(key)) for key in macro_keys])
179+
macro_tensor = np.array([float(eval(key)) for key in macro_keys.values()])
148180

149181
## microscopic body data (particle size, position, etc)
150182
micro_bodies_data = []
151-
for bodyKey in micro_keys_bodies:
183+
for bodyKey in micro_keys_bodies.values():
152184
micro_bodies_data.append([float(eval('b.'+bodyKey)) for b in O.bodies])
153185
micro_bodies_data = np.array(micro_bodies_data)
154186

155187
## macroscopic interaction data (contact pairs, interaction forces, etc)
156188
micro_inters_data = []
157-
if not STORE_EDGE_FEATURES: keys = micro_keys_inters[:2]
189+
if not STORE_EDGE_FEATURES: keys = list(micro_keys_inters.values())[:2]
158190
for key in keys:
159191
micro_inters_data.append([float(eval('i.'+key)) for i in O.interactions if i.isReal])
160192
micro_inters_data = np.array(micro_inters_data)
@@ -167,18 +199,19 @@ def main(pressure, experiment_type,numParticles):
167199
e = np.transpose(e, (1, 0))
168200
n = np.transpose(micro_bodies_data, (1, 0))
169201

170-
f_step = f_sample.create_group(f'{step}')
202+
f_step = f_sample_time.create_group(f'{i}')
171203
f_step['sources'] = src
172204
f_step['destinations'] = dst
173205
if STORE_EDGE_FEATURES: f_step['edge_features'] = e
174-
f_step['node_features'] = n
175-
f_step['macro_input_features'] = macro_tensor[:3]
176-
f_step['macro_output_features'] = macro_tensor[4:7]
177-
f_step['other_features'] = macro_tensor[8:]
206+
f_step['node_features'] = n.astype('float32')
207+
f_step['macro_input_features'] = macro_tensor[0:3]
208+
f_step['macro_output_features'] = macro_tensor[3:6]
209+
f_step['time'] = macro_tensor[6:]
178210

179211
print(f'Added data to {h5file_name}')
180212

213+
181214
for pressure in ['0.1e6']:
182215
for experiment_type in ['drained']:
183-
for numParticles in ['15000']:
216+
for numParticles in ['5000']:
184217
main(pressure, experiment_type, numParticles)

0 commit comments

Comments
 (0)