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 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.001
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
185 final_time = 0.06
186 num_steps = 200000
187 # final_time = 0.00002
188 # num_steps = 2
189 num_outputs = 40
190 dt_out_n = num_steps / num_outputs
191 perform_out = True
192
193
194 poisson1 = 0.25
195 rho1 = 1200.
196 K1 = 2.16e+7
197 E1 = get_E(K1, poisson1)
198 G1 = get_G(E1, poisson1)
199 Gc1 = 50.
200
201 poisson2 = 0.25
202 rho2 = 1200.
203 K2 = 2.16e+7
204 E2 = get_E(K2, poisson2)
205 G2 = get_G(E2, poisson2)
206 Gc2 = 50.
207
208
211 R_contact_factor = 0.95
212
213 # Kn_V_max = 7.385158e+05
214 # Kn = np.power(Kn_V_max, 2)
215 # compute from bulk modulus
216
217 # from bulk modulus
218 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
219 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
220 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
221
222 beta_n_eps = 0.9
223 friction_coeff = 0.5
224 damping_active = False
225 friction_active = False
226 beta_n_factor = 100.
227
228
229 gravity_active = True
230 gravity = [0., -10., 0.]
231
232
233 free_fall_dist = particle_dist - horizon
234 free_fall_vel = [0., 0., 0.]
235 free_fall_vel[1] = -np.sqrt(2. * np.abs(gravity[1]) * free_fall_dist)
236 # free_fall_vel[1] = -2.
237
238
239 neigh_search_factor = 10.
240 neigh_search_interval = 40
241 neigh_search_criteria = "simple_all"
242
243
244
248 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
249 inpf.write("Model:\n")
250 inpf.write(" Dimension: 3\n")
251 inpf.write(" Discretization_Type:\n")
252 inpf.write(" Spatial: finite_difference\n")
253 inpf.write(" Time: central_difference\n")
254 inpf.write(" Final_Time: %4.6e\n" % (final_time))
255 inpf.write(" Time_Steps: %d\n" % (num_steps))
256
257 inpf.write("Policy:\n")
258 inpf.write(" Enable_PostProcessing: true\n")
259
260 #
261 # container info
262 #
263 inpf.write("Container:\n")
264 inpf.write(" Geometry:\n")
265 inpf.write(" Type: cuboid\n")
266 contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + particle_dist, 0.]
267 if R2 > R1:
268 contain_params[3] = 2.*R2
269 inpf.write(" Parameters: " + print_dbl_list(contain_params))
270
271 #
272 # zone info
273 #
274 inpf.write("Zone:\n")
275 inpf.write(" Zones: 2\n")
276
277
278 inpf.write(" Zone_1:\n")
279 inpf.write(" Is_Wall: false\n")
280
281
282 inpf.write(" Zone_2:\n")
283 inpf.write(" Is_Wall: false\n")
284
285 #
286 # particle info
287 #
288 inpf.write("Particle:\n")
289 inpf.write(" Test_Name: two_particle_3D\n")
290 inpf.write(" Zone_1:\n")
291 inpf.write(" Type: sphere\n")
292 p1_geom = [R1, center[0], center[1], center[2]]
293 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
294 inpf.write(" Zone_2:\n")
295 inpf.write(" Type: sphere\n")
296 p2_geom = [R2, center[0], center[1], center[2]]
297 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
298
299 #
300 # particle generation
301 #
302 inpf.write("Particle_Generation:\n")
303 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
304 inpf.write(" File_Data_Type: loc_rad_orient\n")
305
306 #
307 # Mesh info
308 #
309 inpf.write("Mesh:\n")
310
311
312 inpf.write(" Zone_1:\n")
313 inpf.write(" File: mesh_sphere_1_" + str(pp_tag) + ".msh \n")
314
315
316 inpf.write(" Zone_2:\n")
317 inpf.write(" File: mesh_sphere_2_" + str(pp_tag) + ".msh \n")
318
319 #
320 # Contact info
321 #
322 inpf.write("Contact:\n")
323
324
325 inpf.write(" Zone_11:\n")
326 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
327
328 if damping_active == False:
329 inpf.write(" Damping_On: false\n")
330 if friction_active == False:
331 inpf.write(" Friction_On: false\n")
332
333 inpf.write(" Kn: %4.6e\n" % (Kn_11))
334 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
335 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
336 inpf.write(" Kn_Factor: 1.0\n")
337 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
338
339
340 inpf.write(" Zone_12:\n")
341 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
342
343 if damping_active == False:
344 inpf.write(" Damping_On: false\n")
345 if friction_active == False:
346 inpf.write(" Friction_On: false\n")
347
348 inpf.write(" Kn: %4.6e\n" % (Kn_12))
349 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
350 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
351 inpf.write(" Kn_Factor: 1.0\n")
352 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
353
354
355 inpf.write(" Zone_22:\n")
356 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
357
358 if damping_active == False:
359 inpf.write(" Damping_On: false\n")
360 if friction_active == False:
361 inpf.write(" Friction_On: false\n")
362
363 inpf.write(" Kn: %4.6e\n" % (Kn_22))
364 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
365 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
366 inpf.write(" Kn_Factor: 1.0\n")
367 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
368
369 # Neighbor info
370 inpf.write("Neighbor:\n")
371 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
372 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
373 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
374
375 #
376 # Material info
377 #
378 inpf.write("Material:\n")
379
380
381 inpf.write(" Zone_1:\n")
382 inpf.write(" Type: PDState\n")
383 inpf.write(" Horizon: %4.6e\n" % (horizon))
384 inpf.write(" Density: %4.6e\n" % (rho1))
385 inpf.write(" Compute_From_Classical: true\n")
386 inpf.write(" K: %4.6e\n" % (K1))
387 inpf.write(" G: %4.6e\n" % (G1))
388 inpf.write(" Gc: %4.6e\n" % (Gc1))
389 inpf.write(" Influence_Function:\n")
390 inpf.write(" Type: 1\n")
391
392
393 inpf.write(" Zone_2:\n")
394 inpf.write(" Type: PDState\n")
395 inpf.write(" Horizon: %4.6e\n" % (horizon))
396 inpf.write(" Density: %4.6e\n" % (rho2))
397 inpf.write(" Compute_From_Classical: true\n")
398 inpf.write(" K: %4.6e\n" % (K2))
399 inpf.write(" G: %4.6e\n" % (G2))
400 inpf.write(" Gc: %4.6e\n" % (Gc2))
401 inpf.write(" Influence_Function:\n")
402 inpf.write(" Type: 1\n")
403
404 #
405 # Force
406 #
407 if gravity_active == True:
408 inpf.write("Force_BC:\n")
409 inpf.write(" Gravity: " + print_dbl_list(gravity))
410
411 #
412 # IC
413 #
414 inpf.write("IC:\n")
415 inpf.write(" Constant_Velocity:\n")
416 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
417 inpf.write(" Particle_List: [1]\n")
418
419 #
420 # Displacement
421 #
422 inpf.write("Displacement_BC:\n")
423 inpf.write(" Sets: 1\n")
424
425 inpf.write(" Set_1:\n")
426 inpf.write(" Particle_List: [0]\n")
427 inpf.write(" Direction: [1,2,3]\n")
428 inpf.write(" Time_Function:\n")
429 inpf.write(" Type: constant\n")
430 inpf.write(" Parameters:\n")
431 inpf.write(" - 0.0\n")
432 inpf.write(" Spatial_Function:\n")
433 inpf.write(" Type: constant\n")
434 inpf.write(" Zero_Displacement: true\n")
435
436 #
437 # Output info
438 #
439 inpf.write("Output:\n")
440 inpf.write(" Path: ../out/\n")
441 inpf.write(" Tags:\n")
442 inpf.write(" - Displacement\n")
443 inpf.write(" - Velocity\n")
444 inpf.write(" - Force\n")
445 inpf.write(" - Force_Density\n")
446 inpf.write(" - Damage_Z\n")
447 inpf.write(" - Damage\n")
448 inpf.write(" - Nodal_Volume\n")
449 inpf.write(" - Zone_ID\n")
450 inpf.write(" - Particle_ID\n")
451 inpf.write(" - Fixity\n")
452 inpf.write(" - Force_Fixity\n")
453 inpf.write(" - Contact_Nodes\n")
454 inpf.write(" - No_Fail_Node\n")
455 inpf.write(" - Boundary_Node_Flag\n")
456 inpf.write(" - Theta\n")
457 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
458 inpf.write(" Compress_Type: zlib\n")
459 inpf.write(" Perform_FE_Out: false\n")
460 if perform_out:
461 inpf.write(" Perform_Out: true\n")
462 else:
463 inpf.write(" Perform_Out: false\n")
464 inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))
465
466 inpf.write(" Debug: 3\n")
467 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
468
469 # close file
470 inpf.close()
471
472
473 # generate particle locations
474 particle_locations(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
475
476 p_mesh_fname = ['mesh_sphere_1', 'mesh_sphere_2']
477 # generate particle .geo file (large)
478 generate_particle_gmsh_input(inp_dir, 'mesh_sphere_1', center, R1, mesh_size, pp_tag)
479 generate_particle_gmsh_input(inp_dir, 'mesh_sphere_2', center, R2, mesh_size, pp_tag)
480
481 os.system("mkdir -p ../out")
482
483 for s in p_mesh_fname:
484 print('\n\n')
485 print(s)
486 print("gmsh {}_{}.geo -3".format(s, pp_tag))
487 print('\n\n')
488 os.system("gmsh {}_{}.geo -3".format(s, pp_tag))
489 os.system("gmsh {}_{}.geo -3 -o {}_{}.vtk".format(s, pp_tag, s, pp_tag))
490
491
492
494inp_dir = './'
495pp_tag = 0
496if len(sys.argv) > 1:
497 pp_tag = int(sys.argv[1])
498
499create_input_file(inp_dir, pp_tag)