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, R1, R1])
74 sim_particles.append([1., R1, 2. * R1 + R2 + offset, R1, R2])
75
76 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
77 inpf.write("i, x, y, z, r\n")
78 for p in sim_particles:
79 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
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 geof.write("Point(6) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz + sim_radius, sim_h))
113 geof.write("Point(7) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz - sim_radius, sim_h))
114
115 #
116 # circlular arc
117 #
118 geof.write("Circle(1) = {2, 1, 4};\n")
119 geof.write("Circle(2) = {4, 1, 3};\n")
120 geof.write("Circle(3) = {3, 1, 5};\n")
121 geof.write("Circle(4) = {5, 1, 2};\n")
122
123 geof.write("Circle(5) = {2, 1, 6};\n")
124 geof.write("Circle(6) = {6, 1, 3};\n")
125 geof.write("Circle(7) = {3, 1, 7};\n")
126 geof.write("Circle(8) = {7, 1, 2};\n")
127
128 geof.write("Circle(9) = {4, 1, 6};\n")
129 geof.write("Circle(10) = {6, 1, 5};\n")
130 geof.write("Circle(11) = {5, 1, 7};\n")
131 geof.write("Circle(12) = {7, 1, 4};\n")
132
133 #
134 # surfaces
135 #
136 geof.write("Line Loop(14) = {2, 7, 12};\n")
137 geof.write("Ruled Surface(14) = {14};\n")
138
139 geof.write("Line Loop(16) = {2, -6, -9};\n")
140 geof.write("Ruled Surface(16) = {16};\n")
141 geof.write("Line Loop(18) = {3, -10, 6};\n")
142 geof.write("Ruled Surface(18) = {18};\n")
143 geof.write("Line Loop(20) = {3, 11, -7};\n")
144 geof.write("Ruled Surface(20) = {20};\n")
145 geof.write("Line Loop(22) = {4, -8, -11};\n")
146 geof.write("Ruled Surface(22) = {22};\n")
147 geof.write("Line Loop(24) = {4, 5, 10};\n")
148 geof.write("Ruled Surface(24) = {24};\n")
149 geof.write("Line Loop(26) = {1, 9, -5};\n")
150 geof.write("Ruled Surface(26) = {26};\n")
151 geof.write("Line Loop(28) = {1, -12, 8};\n")
152 geof.write("Ruled Surface(28) = {28};\n")
153
154 geof.write("Surface Loop(30) = {14, 16, 18, 20, 22, 24, 26, 28};\n")
155 # tag = '"' + "a" + '"'
156 # geof.write("Physical Surface(%s) = {30};\n" % (tag))
157 geof.write("Volume(30) = {30};\n")
158 geof.write("Physical Volume(1) = {30};\n")
159
160 # add center point to plane surface
161 geof.write("Point{1} In Volume {1};")
162
163 # close file
164 geof.close()
165
166
167
168def create_input_file(inp_dir, pp_tag):
169 """Generates input file for two-particle test"""
170
171 sim_inp_dir = str(inp_dir)
172
173
174 center = [0., 0., 0.]
175 R1 = 0.003
176 R2 = 0.001
177 mesh_size = R1 / 5.
178 if R2 < R1:
179 mesh_size = R2 / 5.
180
181 horizon = 3. * mesh_size
182 particle_dist = 0.001
183
184 # maximum distance to stop the simulation
185 max_dist_check = 20.*R1
186
187
188 final_time = 0.04
189 num_steps = 20000
190 # final_time = 0.00002
191 # num_steps = 2
192 num_outputs = 10
193 dt_out_n = num_steps / num_outputs
194 test_dt_out_n = dt_out_n / 100
195 perform_out = True
196
197
198 poisson1 = 0.25
199 rho1 = 100.
200 K1 = 1e+4
201 E1 = get_E(K1, poisson1)
202 G1 = get_G(E1, poisson1)
203 Gc1 = 10.
204
205 poisson2 = 0.25
206 rho2 = 100.
207 K2 = 1e+4
208 E2 = get_E(K2, poisson2)
209 G2 = get_G(E2, poisson2)
210 Gc2 = 10.
211
212
215 R_contact_factor = 0.95
216
217 # Kn_V_max = 7.385158e+05
218 # Kn = np.power(Kn_V_max, 2)
219 # compute from bulk modulus
220
221 # from bulk modulus (power 6 for dimension 3)
222 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 6))
223 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 6))
224 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 6))
225
226 beta_n_eps = 0.95
227 friction_coeff = 0.5
228 damping_active = False
229 friction_active = False
230 beta_n_factor = 10.
231
232
233 gravity_active = True
234 gravity = [0., -10., 0.]
235
236
237 free_fall_dist = particle_dist - horizon
238 free_fall_vel = [0., 0., 0.]
239 free_fall_vel[1] = -np.sqrt(2. * np.abs(gravity[1]) * free_fall_dist)
240 #free_fall_vel[1] = -1.
241
242
243
247 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
248 inpf.write("Model:\n")
249 inpf.write(" Dimension: 3\n")
250 inpf.write(" Discretization_Type:\n")
251 inpf.write(" Spatial: finite_difference\n")
252 inpf.write(" Time: central_difference\n")
253 inpf.write(" Final_Time: %4.6e\n" % (final_time))
254 inpf.write(" Time_Steps: %d\n" % (num_steps))
255
256 inpf.write("Policy:\n")
257 inpf.write(" Enable_PostProcessing: true\n")
258
259 #
260 # zone info
261 #
262 inpf.write("Zone:\n")
263 inpf.write(" Zones: 2\n")
264
265
266 inpf.write(" Zone_1:\n")
267 inpf.write(" Is_Wall: false\n")
268 inpf.write(" Type: sphere\n")
269 zone_1_circle = [R1, center[0], center[1], center[2]]
270 inpf.write(" Parameters: " + print_dbl_list(zone_1_circle))
271
272
273 inpf.write(" Zone_2:\n")
274 inpf.write(" Is_Wall: false\n")
275 inpf.write(" Type: sphere\n")
276 zone_2_circle = [R2, center[0], center[1], center[2]]
277 inpf.write(" Parameters: " + print_dbl_list(zone_2_circle))
278
279 #
280 # particle info
281 #
282 inpf.write("Particle:\n")
283 inpf.write(" Test_Name: two_particle\n")
284 inpf.write(" Zone_1:\n")
285 inpf.write(" Type: sphere\n")
286 inpf.write(" Parameters: [%4.6e]\n" % (R1))
287 inpf.write(" Zone_2:\n")
288 inpf.write(" Type: sphere\n")
289 inpf.write(" Parameters: [%4.6e]\n" % (R2))
290
291 #
292 # particle generation
293 #
294 inpf.write("Particle_Generation:\n")
295 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
296
297 #
298 # container info
299 #
300 inpf.write("Container:\n")
301 inpf.write(" Geometry:\n")
302 inpf.write(" Type: cuboid\n")
303 contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + particle_dist, 2.*R1]
304 if R2 > R1:
305 contain_params[3] = 2.*R2
306 contain_params[5] = 2.*R2
307 inpf.write(" Parameters: " + print_dbl_list(contain_params))
308
309 #
310 # Mesh info
311 #
312 inpf.write("Mesh:\n")
313
314
315 inpf.write(" Zone_1:\n")
316 inpf.write(" File: mesh_cir_1_" + str(pp_tag) + ".msh \n")
317 inpf.write(" Reference_Particle:\n")
318 inpf.write(" Type: sphere\n")
319 inpf.write(" Parameters: [%4.6e]\n" % (R1))
320
321
322 inpf.write(" Zone_2:\n")
323 inpf.write(" File: mesh_cir_2_" + str(pp_tag) + ".msh \n")
324 inpf.write(" Reference_Particle:\n")
325 inpf.write(" Type: sphere\n")
326 inpf.write(" Parameters: [%4.6e]\n" % (R2))
327
328 #
329 # Contact info
330 #
331 inpf.write("Contact:\n")
332
333
334 inpf.write(" Zone_11:\n")
335 # inpf.write(" Contact_Radius: %4.6e\n" % (R_contact))
336 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
337
338 if damping_active == False:
339 inpf.write(" Damping_On: false\n")
340 if friction_active == False:
341 inpf.write(" Friction_On: false\n")
342
343 inpf.write(" Kn: %4.6e\n" % (Kn_11))
344 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
345 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
346 inpf.write(" Kn_Factor: 1.0\n")
347 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
348
349
350 inpf.write(" Zone_12:\n")
351 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
352
353 if damping_active == False:
354 inpf.write(" Damping_On: false\n")
355 if friction_active == False:
356 inpf.write(" Friction_On: false\n")
357
358 inpf.write(" Kn: %4.6e\n" % (Kn_12))
359 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
360 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
361 inpf.write(" Kn_Factor: 1.0\n")
362 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
363
364
365 inpf.write(" Zone_22:\n")
366 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
367
368 if damping_active == False:
369 inpf.write(" Damping_On: false\n")
370 if friction_active == False:
371 inpf.write(" Friction_On: false\n")
372
373 inpf.write(" Kn: %4.6e\n" % (Kn_22))
374 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
375 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
376 inpf.write(" Kn_Factor: 1.0\n")
377 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
378
379 #
380 # Neighbor info
381 #
382 inpf.write("Neighbor:\n")
383 inpf.write(" Update_Criteria: simple_all\n")
384 inpf.write(" Search_Factor: 5.0\n")
385
386 #
387 # Material info
388 #
389 inpf.write("Material:\n")
390
391
392 inpf.write(" Zone_1:\n")
393 inpf.write(" Type: PDState\n")
394 inpf.write(" Horizon: %4.6e\n" % (horizon))
395 inpf.write(" Density: %4.6e\n" % (rho1))
396 inpf.write(" Compute_From_Classical: true\n")
397 inpf.write(" K: %4.6e\n" % (K1))
398 inpf.write(" G: %4.6e\n" % (G1))
399 inpf.write(" Gc: %4.6e\n" % (Gc1))
400 inpf.write(" Influence_Function:\n")
401 inpf.write(" Type: 1\n")
402
403
404 inpf.write(" Zone_2:\n")
405 inpf.write(" Type: PDState\n")
406 inpf.write(" Horizon: %4.6e\n" % (horizon))
407 inpf.write(" Density: %4.6e\n" % (rho2))
408 inpf.write(" Compute_From_Classical: true\n")
409 inpf.write(" K: %4.6e\n" % (K2))
410 inpf.write(" G: %4.6e\n" % (G2))
411 inpf.write(" Gc: %4.6e\n" % (Gc2))
412 inpf.write(" Influence_Function:\n")
413 inpf.write(" Type: 1\n")
414
415 #
416 # Force
417 #
418 if gravity_active == True:
419 inpf.write("Force_BC:\n")
420 inpf.write(" Gravity: " + print_dbl_list(gravity))
421
422 #
423 # IC
424 #
425 inpf.write("IC:\n")
426 inpf.write(" Constant_Velocity:\n")
427 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
428 inpf.write(" Particle_List: [1]\n")
429
430 #
431 # Displacement
432 #
433 inpf.write("Displacement_BC:\n")
434 inpf.write(" Sets: 1\n")
435
436 inpf.write(" Set_1:\n")
437 inpf.write(" Particle_List: [0]\n")
438 inpf.write(" Direction: [1,2,3]\n")
439 inpf.write(" Time_Function:\n")
440 inpf.write(" Type: constant\n")
441 inpf.write(" Parameters:\n")
442 inpf.write(" - 0.0\n")
443 inpf.write(" Spatial_Function:\n")
444 inpf.write(" Type: constant\n")
445 inpf.write(" Zero_Displacement: true\n")
446
447 #
448 # Output info
449 #
450 inpf.write("Output:\n")
451 inpf.write(" Path: ../out/\n")
452 inpf.write(" Tags:\n")
453 inpf.write(" - Displacement\n")
454 inpf.write(" - Velocity\n")
455 inpf.write(" - Force\n")
456 inpf.write(" - Force_Density\n")
457 inpf.write(" - Damage_Z\n")
458 inpf.write(" - Damage\n")
459 inpf.write(" - Nodal_Volume\n")
460 inpf.write(" - Zone_ID\n")
461 inpf.write(" - Particle_ID\n")
462 inpf.write(" - Fixity\n")
463 inpf.write(" - Force_Fixity\n")
464 inpf.write(" - Contact_Nodes\n")
465 inpf.write(" - No_Fail_Node\n")
466 inpf.write(" - Boundary_Node_Flag\n")
467 inpf.write(" - Theta\n")
468 inpf.write(" - Contact_Data\n")
469 # inpf.write(" - Strain_Stress\n")
470 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
471 inpf.write(" Compress_Type: zlib\n")
472 inpf.write(" Perform_FE_Out: false\n")
473 if perform_out:
474 inpf.write(" Perform_Out: true\n")
475 else:
476 inpf.write(" Perform_Out: false\n")
477 inpf.write(" Test_Output_Interval: %d\n" % (test_dt_out_n))
478
479 inpf.write(" Debug: 1\n")
480 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
481 inpf.write(" Output_Criteria: \n")
482 # inpf.write(" Type: max_particle_dist\n")
483 # inpf.write(" Parameters: [%4.6e]\n" % (2. * sim_h))
484 inpf.write(" Type: max_node_dist\n")
485 inpf.write(" Parameters: [%4.6e]\n" % (max_dist_check))
486
487
488 inpf.write("HPX:\n")
489 inpf.write(" Partitions: 1\n")
490
491 # close file
492 inpf.close()
493
494
495 # generate particle locations
496 particle_locations(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
497
498 # generate particle .geo file (large)
499 generate_particle_gmsh_input(inp_dir, 'mesh_cir_1', center, R1, mesh_size, pp_tag)
500 generate_particle_gmsh_input(inp_dir, 'mesh_cir_2', center, R2, mesh_size, pp_tag)
501
502
503
505inp_dir = './'
506pp_tag = 0
507if len(sys.argv) > 1:
508 pp_tag = int(sys.argv[1])
509
510create_input_file(inp_dir, pp_tag)