PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
problem_setup.py
Go to the documentation of this file.
1import os
2import numpy as np
3import sys
4
5def print_bool(arg, prefix = ""):
6
7 str = prefix
8 if arg == True:
9 str += "True\n"
10 else:
11 str += "False\n"
12 return str
13
14def print_dbl(arg, prefix = ""):
15
16 str = prefix + "%4.6e\n" % (arg)
17 return str
18
19def print_int(arg, prefix = ""):
20 str = prefix + "%d\n" % (arg)
21 return str
22
23def print_dbl_list(arg, prefix = ""):
24 str = prefix + "["
25 N = len(arg)
26 for i in range(N):
27 str += "%4.6e" % (arg[i])
28 if i < N - 1:
29 str += ", "
30 else:
31 str += "]\n"
32
33 return str
34
35def print_int_list(arg, prefix = ""):
36 str = prefix + "["
37 N = len(arg)
38 for i in range(N):
39 str += "%d" % (arg[i])
40 if i < N - 1:
41 str += ", "
42 else:
43 str += "]\n"
44
45 return str
46
47def does_intersect(p, r, R, particles, padding):
48
49 for q in particles:
50
51 pq = np.array([p[i] - q[i] for i in range(3)])
52 if np.linalg.norm(pq) <= r + R + padding:
53 return True
54
55 return False
56
57
58def get_E(K, nu):
59 return 3. * K * (1. - 2. * nu)
60
61def get_G(E, nu):
62 return E / (2. * (1. + nu))
63
64
65def get_eff_k(k1, k2):
66 return 2. * k1 * k2 / (k1 + k2)
67
68
69def particle_locations(inp_dir, pp_tag, R1, R2, offset):
70 """Generate particle location data"""
71
72 sim_particles = []
73 sim_particles.append([0., R1, R1, 0., R1, 0.])
74 sim_particles.append([1., R1, 2. * R1 + R2 + offset, 0., R2, np.pi*0.5])
75
76 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
77 inpf.write("i, x, y, z, r, o\n")
78 for p in sim_particles:
79 inpf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4], p[5]))
80
81 inpf.close()
82
83
84def generate_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
85
86 sim_inp_dir = str(inp_dir)
87
88 # center and radius
89 sim_Cx = center[0]
90 sim_Cy = center[1]
91 sim_Cz = center[2]
92 sim_radius = radius
93
94 # mesh size
95 sim_h = mesh_size
96
97 #
98 # create .geo file for gmsh
99 #
100 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
101 geof.write("cl__1 = 1;\n")
102 geof.write("Mesh.MshFileVersion = 2.2;\n")
103
104 #
105 # points
106 #
107 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz, sim_h));
108 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx + sim_radius, sim_Cy, sim_Cz, sim_h))
109 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx - sim_radius, sim_Cy, sim_Cz, sim_h))
110 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy + sim_radius, sim_Cz, sim_h))
111 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy - sim_radius, sim_Cz, sim_h))
112
113 #
114 # circlular arc
115 #
116 geof.write("Circle(1) = {2, 1, 4};\n")
117 geof.write("Circle(2) = {4, 1, 3};\n")
118 geof.write("Circle(3) = {3, 1, 5};\n")
119 geof.write("Circle(4) = {5, 1, 2};\n")
120
121 #
122 # surfaces
123 #
124 geof.write("Line Loop(1) = {2, 3, 4, 1};\n")
125
126 #
127 # plane surface
128 #
129 geof.write("Plane Surface(1) = {1};\n")
130
131 #
132 # physical surface
133 #
134 # tag = '"' + "a" + '"'
135 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
136
137 # # add center point to plane surface
138 geof.write("Point{1} In Surface {1};")
139
140 # close file
141 geof.close()
142
143
144
145def create_input_file(inp_dir, pp_tag):
146 """Generates input file for two-particle test"""
147
148 sim_inp_dir = str(inp_dir)
149
150
151 center = [0., 0., 0.]
152 R1 = 0.001
153 R2 = 0.001
154 mesh_size = R1 / 5.
155 if R2 < R1:
156 mesh_size = R2 / 5.
157
158 horizon = 3. * mesh_size
159 particle_dist = 0.001
160
161
162 final_time = 0.04
163 num_steps = 200000
164 # final_time = 0.00002
165 # num_steps = 2
166 num_outputs = 100
167 dt_out_n = num_steps / num_outputs
168 perform_out = True
169
170
171 poisson1 = 0.25
172 rho1 = 1200.
173 K1 = 2.16e+7
174 E1 = get_E(K1, poisson1)
175 G1 = get_G(E1, poisson1)
176 Gc1 = 50.
177
178 poisson2 = 0.25
179 rho2 = 1200.
180 K2 = 2.16e+7
181 E2 = get_E(K2, poisson2)
182 G2 = get_G(E2, poisson2)
183 Gc2 = 50.
184
185
188 R_contact_factor = 0.95
189
190 # Kn_V_max = 7.385158e+05
191 # Kn = np.power(Kn_V_max, 2)
192 # compute from bulk modulus
193
194 # from bulk modulus
195 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
196 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
197 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
198
199 beta_n_eps = 0.9
200 friction_coeff = 0.5
201 damping_active = True
202 friction_active = False
203 beta_n_factor = 100.
204
205
206 gravity_active = True
207 gravity = [0., -10., 0.]
208
209
210 free_fall_dist = particle_dist - horizon
211 free_fall_vel = [0., 0., 0.]
212 free_fall_vel[1] = -np.sqrt(2. * np.abs(gravity[1]) * free_fall_dist)
213 # free_fall_vel[1] = -2.
214
215
216
220 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
221 inpf.write("Model:\n")
222 inpf.write(" Dimension: 2\n")
223 inpf.write(" Discretization_Type:\n")
224 inpf.write(" Spatial: finite_difference\n")
225 inpf.write(" Time: central_difference\n")
226 inpf.write(" Final_Time: %4.6e\n" % (final_time))
227 inpf.write(" Time_Steps: %d\n" % (num_steps))
228
229 inpf.write("Policy:\n")
230 inpf.write(" Enable_PostProcessing: true\n")
231
232 #
233 # container info
234 #
235 inpf.write("Container:\n")
236 inpf.write(" Geometry:\n")
237 inpf.write(" Type: rectangle\n")
238 contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + particle_dist, 0.]
239 if R2 > R1:
240 contain_params[3] = 2.*R2
241 inpf.write(" Parameters: " + print_dbl_list(contain_params))
242
243 #
244 # zone info
245 #
246 inpf.write("Zone:\n")
247 inpf.write(" Zones: 2\n")
248
249
250 inpf.write(" Zone_1:\n")
251 inpf.write(" Is_Wall: false\n")
252
253
254 inpf.write(" Zone_2:\n")
255 inpf.write(" Is_Wall: false\n")
256
257 #
258 # particle info
259 #
260 inpf.write("Particle:\n")
261 inpf.write(" Test_Name: two_particle\n")
262 inpf.write(" Zone_1:\n")
263 inpf.write(" Type: circle\n")
264 p1_geom = [R1, center[0], center[1], center[2]]
265 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
266 inpf.write(" Zone_2:\n")
267 inpf.write(" Type: circle\n")
268 p2_geom = [R2, center[0], center[1], center[2]]
269 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
270
271 #
272 # particle generation
273 #
274 inpf.write("Particle_Generation:\n")
275 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
276 inpf.write(" File_Data_Type: loc_rad_orient\n")
277
278 #
279 # Mesh info
280 #
281 inpf.write("Mesh:\n")
282
283
284 inpf.write(" Zone_1:\n")
285 inpf.write(" File: mesh_cir_1_" + str(pp_tag) + ".msh \n")
286
287
288 inpf.write(" Zone_2:\n")
289 inpf.write(" File: mesh_cir_2_" + str(pp_tag) + ".msh \n")
290
291 #
292 # Contact info
293 #
294 inpf.write("Contact:\n")
295
296
297 inpf.write(" Zone_11:\n")
298 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
299
300 if damping_active == False:
301 inpf.write(" Damping_On: false\n")
302 if friction_active == False:
303 inpf.write(" Friction_On: false\n")
304
305 inpf.write(" Kn: %4.6e\n" % (Kn_11))
306 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
307 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
308 inpf.write(" Kn_Factor: 1.0\n")
309 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
310
311
312 inpf.write(" Zone_12:\n")
313 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
314
315 if damping_active == False:
316 inpf.write(" Damping_On: false\n")
317 if friction_active == False:
318 inpf.write(" Friction_On: false\n")
319
320 inpf.write(" Kn: %4.6e\n" % (Kn_12))
321 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
322 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
323 inpf.write(" Kn_Factor: 1.0\n")
324 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
325
326
327 inpf.write(" Zone_22:\n")
328 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
329
330 if damping_active == False:
331 inpf.write(" Damping_On: false\n")
332 if friction_active == False:
333 inpf.write(" Friction_On: false\n")
334
335 inpf.write(" Kn: %4.6e\n" % (Kn_22))
336 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
337 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
338 inpf.write(" Kn_Factor: 1.0\n")
339 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
340
341 #
342 # Neighbor info
343 #
344 inpf.write("Neighbor:\n")
345 inpf.write(" Update_Criteria: simple_all\n")
346 inpf.write(" Search_Factor: 5.0\n")
347
348 #
349 # Material info
350 #
351 inpf.write("Material:\n")
352
353
354 inpf.write(" Zone_1:\n")
355 inpf.write(" Type: PDState\n")
356 inpf.write(" Horizon: %4.6e\n" % (horizon))
357 inpf.write(" Density: %4.6e\n" % (rho1))
358 inpf.write(" Compute_From_Classical: true\n")
359 inpf.write(" K: %4.6e\n" % (K1))
360 inpf.write(" G: %4.6e\n" % (G1))
361 inpf.write(" Gc: %4.6e\n" % (Gc1))
362 inpf.write(" Influence_Function:\n")
363 inpf.write(" Type: 1\n")
364
365
366 inpf.write(" Zone_2:\n")
367 inpf.write(" Type: PDState\n")
368 inpf.write(" Horizon: %4.6e\n" % (horizon))
369 inpf.write(" Density: %4.6e\n" % (rho2))
370 inpf.write(" Compute_From_Classical: true\n")
371 inpf.write(" K: %4.6e\n" % (K2))
372 inpf.write(" G: %4.6e\n" % (G2))
373 inpf.write(" Gc: %4.6e\n" % (Gc2))
374 inpf.write(" Influence_Function:\n")
375 inpf.write(" Type: 1\n")
376
377 #
378 # Force
379 #
380 if gravity_active == True:
381 inpf.write("Force_BC:\n")
382 inpf.write(" Gravity: " + print_dbl_list(gravity))
383
384 #
385 # IC
386 #
387 inpf.write("IC:\n")
388 inpf.write(" Constant_Velocity:\n")
389 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
390 inpf.write(" Particle_List: [1]\n")
391
392 #
393 # Displacement
394 #
395 inpf.write("Displacement_BC:\n")
396 inpf.write(" Sets: 1\n")
397
398 inpf.write(" Set_1:\n")
399 inpf.write(" Particle_List: [0]\n")
400 inpf.write(" Direction: [1,2]\n")
401 inpf.write(" Time_Function:\n")
402 inpf.write(" Type: constant\n")
403 inpf.write(" Parameters:\n")
404 inpf.write(" - 0.0\n")
405 inpf.write(" Spatial_Function:\n")
406 inpf.write(" Type: constant\n")
407 inpf.write(" Zero_Displacement: true\n")
408
409 #
410 # Output info
411 #
412 inpf.write("Output:\n")
413 inpf.write(" Path: ../out/\n")
414 inpf.write(" Tags:\n")
415 inpf.write(" - Displacement\n")
416 inpf.write(" - Velocity\n")
417 inpf.write(" - Force\n")
418 inpf.write(" - Force_Density\n")
419 inpf.write(" - Damage_Z\n")
420 inpf.write(" - Damage\n")
421 inpf.write(" - Nodal_Volume\n")
422 inpf.write(" - Zone_ID\n")
423 inpf.write(" - Particle_ID\n")
424 inpf.write(" - Fixity\n")
425 inpf.write(" - Force_Fixity\n")
426 inpf.write(" - Contact_Nodes\n")
427 inpf.write(" - No_Fail_Node\n")
428 inpf.write(" - Boundary_Node_Flag\n")
429 inpf.write(" - Theta\n")
430 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
431 inpf.write(" Compress_Type: zlib\n")
432 inpf.write(" Perform_FE_Out: false\n")
433 if perform_out:
434 inpf.write(" Perform_Out: true\n")
435 else:
436 inpf.write(" Perform_Out: false\n")
437 inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))
438
439 inpf.write(" Debug: 1\n")
440 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
441
442
443 inpf.write("HPX:\n")
444 inpf.write(" Partitions: 1\n")
445
446 # close file
447 inpf.close()
448
449
450 # generate particle locations
451 particle_locations(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
452
453 # generate particle .geo file (large)
454 generate_particle_gmsh_input(inp_dir, 'mesh_cir_1', center, R1, mesh_size, pp_tag)
455 generate_particle_gmsh_input(inp_dir, 'mesh_cir_2', center, R2, mesh_size, pp_tag)
456
457
458
460inp_dir = './'
461pp_tag = 0
462if len(sys.argv) > 1:
463 pp_tag = int(sys.argv[1])
464
465create_input_file(inp_dir, pp_tag)
problem_setup.get_E
get_E(K, nu)
Definition problem_setup.py:136
problem_setup.get_eff_k
get_eff_k(k1, k2)
Definition problem_setup.py:143
problem_setup.generate_particle_gmsh_input
generate_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag)
Definition problem_setup.py:84
problem_setup.create_input_file
create_input_file(inp_dir, pp_tag)
Definition problem_setup.py:793
problem_setup.particle_locations
particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R2, id_choices1, id_choices2, Nmax, R_in, bar_rect, z_coord, add_orient=True)
Definition problem_setup.py:287
problem_setup.print_dbl_list
print_dbl_list(arg, prefix="")
Definition problem_setup.py:25
problem_setup.get_G
get_G(E, nu)
Definition problem_setup.py:139