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, wall_particle_dist, wall1_cuboid):
70 """Generate particle location data"""
71
72 sim_particles = []
73 sim_particles.append([0., 0., R1+wall_particle_dist, 0., R1, 0.])
74 sim_particles.append([1., 0.5*(wall1_cuboid[0] + wall1_cuboid[3]), 0.5*(wall1_cuboid[1] + wall1_cuboid[4]), 0.5*(wall1_cuboid[2] + wall1_cuboid[5]), wall1_cuboid[3] - wall1_cuboid[0], 0.])
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_sphere_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
166def generate_cuboid_wall_gmsh_input(inp_dir, filename, cuboid, mesh_size, pp_tag):
167
168 sim_inp_dir = str(inp_dir)
169
170 # center and radius
171 sim_Cx = 0.5*(cuboid[0] + cuboid[3])
172 sim_Cy = 0.5*(cuboid[1] + cuboid[4])
173 sim_Cz = 0.5*(cuboid[2] + cuboid[5])
174 sim_Lx, sim_Ly, sim_Lz = cuboid[3] - cuboid[0], cuboid[4] - cuboid[1], cuboid[5] - cuboid[2]
175
176 sim_corner_p1 = [cuboid[0], cuboid[1], cuboid[2]]
177 sim_corner_p2 = [cuboid[3], cuboid[1], cuboid[2]]
178
179 # mesh size
180 sim_h = mesh_size
181
182 #
183 # create .geo file for gmsh
184 #
185 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
186 geof.write("cl__1 = 1;\n")
187 geof.write("Mesh.MshFileVersion = 2.2;\n")
188
189 #
190 # points
191 #
192 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_corner_p1[0], sim_corner_p1[1], sim_corner_p1[2], sim_h))
193 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_corner_p1[0]+sim_Lx, sim_corner_p1[1], sim_corner_p1[2], sim_h))
194 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_corner_p1[0]+sim_Lx, sim_corner_p1[1]+sim_Ly, sim_corner_p1[2], sim_h))
195 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_corner_p1[0], sim_corner_p1[1]+sim_Ly, sim_corner_p1[2], sim_h))
196
197
198 #
199 # Line
200 #
201 geof.write("Line(1) = {4,3};\n")
202 geof.write("Line(2) = {3,2};\n")
203 geof.write("Line(3) = {2,1};\n")
204 geof.write("Line(4) = {1,4};\n")
205
206 #
207 # surfaces
208 #
209 geof.write("Line Loop(5) = {2,3,4,1};\n")
210 geof.write("Plane Surface(6) = {5};\n")
211 geof.write("Extrude {0, 0.0, %4.6e}{ Surface {6}; }\n" % (sim_Lz))
212
213 # close file
214 geof.close()
215
216
217def create_input_file(inp_dir, pp_tag):
218 """Generates input file for two-particle test"""
219
220 sim_inp_dir = str(inp_dir)
221
222
223 center = [0., 0., 0.]
224 R1 = 0.001
225 mesh_size = R1 / 5.
226
227 horizon = 3. * mesh_size
228 wall_particle_dist = 0.2*R1
229
230 # corner points
231 wall1_cuboid = [-3*R1, -2*horizon, -2*R1, 3*R1, 0, 2*R1]
232 wall1_cuboid_ref = [-3*R1, -horizon, -2*R1, 3*R1, horizon, 2*R1]
233
234
235 final_time = 0.0008
236 num_steps = 100000
237 # final_time = 0.00002
238 # num_steps = 2
239 num_outputs = 100
240 dt_out_n = num_steps / num_outputs
241 perform_out = True
242
243
244 poisson1 = 0.25
245 rho1 = 1200.
246 K1 = 2.16e+9
247 E1 = get_E(K1, poisson1)
248 G1 = get_G(E1, poisson1)
249 Gc1 = 50.
250
251 poisson2 = 0.25
252 rho2 = 1200.
253 K2 = 2.16e+7
254 E2 = get_E(K2, poisson2)
255 G2 = get_G(E2, poisson2)
256 Gc2 = 50.
257
258
261 R_contact_factor = 0.95
262
263 # Kn_V_max = 7.385158e+05
264 # Kn = np.power(Kn_V_max, 2)
265 # compute from bulk modulus
266
267 # from bulk modulus
268 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
269 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
270 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
271
272 beta_n_eps = 0.9
273 friction_coeff = 0.5
274 damping_active = False
275 friction_active = False
276 beta_n_factor = 100.
277
278
279 gravity_active = True
280 gravity = [0., -10., 0.]
281
282
283 free_fall_vel = [0., -2., 0.]
284
285
286 neigh_search_factor = 10.
287 neigh_search_interval = 40
288 neigh_search_criteria = "simple_all"
289
290
291
295 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
296 inpf.write("Model:\n")
297 inpf.write(" Dimension: 3\n")
298 inpf.write(" Discretization_Type:\n")
299 inpf.write(" Spatial: finite_difference\n")
300 inpf.write(" Time: central_difference\n")
301 inpf.write(" Final_Time: %4.6e\n" % (final_time))
302 inpf.write(" Time_Steps: %d\n" % (num_steps))
303
304 inpf.write("Policy:\n")
305 inpf.write(" Enable_PostProcessing: true\n")
306
307 #
308 # container info
309 #
310 inpf.write("Container:\n")
311 inpf.write(" Geometry:\n")
312 inpf.write(" Type: cuboid\n")
313 contain_params = [-2*R1, -2*horizon, -2*R1, 2*R1, 3*R1, 2*R1]
314 inpf.write(" Parameters: " + print_dbl_list(contain_params))
315
316 #
317 # zone info
318 #
319 inpf.write("Zone:\n")
320 inpf.write(" Zones: 2\n")
321
322
323 inpf.write(" Zone_1:\n")
324 inpf.write(" Is_Wall: false\n")
325
326
327 inpf.write(" Zone_2:\n")
328 inpf.write(" Is_Wall: false\n")
329
330 #
331 # particle info
332 #
333 inpf.write("Particle:\n")
334 inpf.write(" Test_Name: sphere_wall_impact\n")
335 inpf.write(" Zone_1:\n")
336 inpf.write(" Type: sphere\n")
337 p1_geom = [R1, center[0], center[1], center[2]]
338 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
339 inpf.write(" Zone_2:\n")
340 inpf.write(" Type: cuboid\n")
341 inpf.write(" Parameters: " + print_dbl_list(wall1_cuboid_ref))
342
343 #
344 # particle generation
345 #
346 inpf.write("Particle_Generation:\n")
347 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
348 inpf.write(" File_Data_Type: loc_rad_orient\n")
349
350 #
351 # Mesh info
352 #
353 inpf.write("Mesh:\n")
354
355
356 inpf.write(" Zone_1:\n")
357 inpf.write(" File: mesh_sphere_" + str(pp_tag) + ".msh \n")
358
359
360 inpf.write(" Zone_2:\n")
361 inpf.write(" File: mesh_cuboid_" + str(pp_tag) + ".msh \n")
362
363 #
364 # Contact info
365 #
366 inpf.write("Contact:\n")
367
368
369 inpf.write(" Zone_11:\n")
370 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
371
372 if damping_active == False:
373 inpf.write(" Damping_On: false\n")
374 if friction_active == False:
375 inpf.write(" Friction_On: false\n")
376
377 inpf.write(" Kn: %4.6e\n" % (Kn_11))
378 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
379 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
380 inpf.write(" Kn_Factor: 1.0\n")
381 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
382
383
384 inpf.write(" Zone_12:\n")
385 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
386
387 if damping_active == False:
388 inpf.write(" Damping_On: false\n")
389 if friction_active == False:
390 inpf.write(" Friction_On: false\n")
391
392 inpf.write(" Kn: %4.6e\n" % (Kn_12))
393 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
394 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
395 inpf.write(" Kn_Factor: 1.0\n")
396 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
397
398
399 inpf.write(" Zone_22:\n")
400 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
401
402 if damping_active == False:
403 inpf.write(" Damping_On: false\n")
404 if friction_active == False:
405 inpf.write(" Friction_On: false\n")
406
407 inpf.write(" Kn: %4.6e\n" % (Kn_22))
408 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
409 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
410 inpf.write(" Kn_Factor: 1.0\n")
411 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
412
413 # Neighbor info
414 inpf.write("Neighbor:\n")
415 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
416 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
417 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
418
419 #
420 # Material info
421 #
422 inpf.write("Material:\n")
423
424
425 inpf.write(" Zone_1:\n")
426 inpf.write(" Type: PDState\n")
427 inpf.write(" Horizon: %4.6e\n" % (horizon))
428 inpf.write(" Density: %4.6e\n" % (rho1))
429 inpf.write(" Compute_From_Classical: true\n")
430 inpf.write(" K: %4.6e\n" % (K1))
431 inpf.write(" G: %4.6e\n" % (G1))
432 inpf.write(" Gc: %4.6e\n" % (Gc1))
433 inpf.write(" Influence_Function:\n")
434 inpf.write(" Type: 1\n")
435
436
437 inpf.write(" Zone_2:\n")
438 inpf.write(" Type: PDState\n")
439 inpf.write(" Horizon: %4.6e\n" % (horizon))
440 inpf.write(" Density: %4.6e\n" % (rho2))
441 inpf.write(" Compute_From_Classical: true\n")
442 inpf.write(" K: %4.6e\n" % (K2))
443 inpf.write(" G: %4.6e\n" % (G2))
444 inpf.write(" Gc: %4.6e\n" % (Gc2))
445 inpf.write(" Influence_Function:\n")
446 inpf.write(" Type: 1\n")
447
448 #
449 # Force
450 #
451 if gravity_active == True:
452 inpf.write("Force_BC:\n")
453 inpf.write(" Gravity: " + print_dbl_list(gravity))
454
455 #
456 # IC
457 #
458 inpf.write("IC:\n")
459 inpf.write(" Constant_Velocity:\n")
460 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
461 inpf.write(" Particle_List: [0]\n")
462
463 #
464 # Displacement
465 #
466 # clamp left and right parts
467 wall1_left = [-2*R1, -2*horizon, -2*R1, -2*R1+horizon, 0., 2*R1]
468 wall1_right = [2*R1-horizon, -2*horizon, -2*R1, 2*R1, 0., 2*R1]
469
470 inpf.write("Displacement_BC:\n")
471 inpf.write(" Sets: 2\n")
472
473 inpf.write(" Set_1:\n")
474 inpf.write(" Region:\n")
475 inpf.write(" Geometry:\n")
476 inpf.write(" Type: cuboid\n")
477 inpf.write(" Parameters: " + print_dbl_list(wall1_left))
478 inpf.write(" Particle_List: [1]\n")
479 inpf.write(" Direction: [1,2,3]\n")
480 inpf.write(" Time_Function:\n")
481 inpf.write(" Type: constant\n")
482 inpf.write(" Parameters:\n")
483 inpf.write(" - 0.0\n")
484 inpf.write(" Spatial_Function:\n")
485 inpf.write(" Type: constant\n")
486 inpf.write(" Zero_Displacement: true\n")
487
488 inpf.write(" Set_2:\n")
489 inpf.write(" Region:\n")
490 inpf.write(" Geometry:\n")
491 inpf.write(" Type: cuboid\n")
492 inpf.write(" Parameters: " + print_dbl_list(wall1_right))
493 inpf.write(" Particle_List: [1]\n")
494 inpf.write(" Direction: [1,2,3]\n")
495 inpf.write(" Time_Function:\n")
496 inpf.write(" Type: constant\n")
497 inpf.write(" Parameters:\n")
498 inpf.write(" - 0.0\n")
499 inpf.write(" Spatial_Function:\n")
500 inpf.write(" Type: constant\n")
501 inpf.write(" Zero_Displacement: true\n")
502
503 #
504 # Output info
505 #
506 inpf.write("Output:\n")
507 inpf.write(" Path: ../out/\n")
508 inpf.write(" Tags:\n")
509 inpf.write(" - Displacement\n")
510 inpf.write(" - Velocity\n")
511 inpf.write(" - Force\n")
512 inpf.write(" - Force_Density\n")
513 inpf.write(" - Damage_Z\n")
514 inpf.write(" - Damage\n")
515 inpf.write(" - Nodal_Volume\n")
516 inpf.write(" - Zone_ID\n")
517 inpf.write(" - Particle_ID\n")
518 inpf.write(" - Fixity\n")
519 inpf.write(" - Force_Fixity\n")
520 inpf.write(" - Contact_Nodes\n")
521 inpf.write(" - No_Fail_Node\n")
522 inpf.write(" - Boundary_Node_Flag\n")
523 inpf.write(" - Theta\n")
524 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
525 inpf.write(" Compress_Type: zlib\n")
526 inpf.write(" Perform_FE_Out: true\n")
527 if perform_out:
528 inpf.write(" Perform_Out: true\n")
529 else:
530 inpf.write(" Perform_Out: false\n")
531 inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))
532
533 inpf.write(" Debug: 3\n")
534 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
535
536 # close file
537 inpf.close()
538
539
540 # generate particle locations
541 particle_locations(inp_dir, pp_tag, R1, wall_particle_dist, wall1_cuboid)
542
543 p_mesh_fname = ['mesh_sphere', 'mesh_cuboid']
544 # generate particle .geo file (large)
545 generate_sphere_particle_gmsh_input(inp_dir, 'mesh_sphere', center, R1, mesh_size, pp_tag)
546 generate_cuboid_wall_gmsh_input(inp_dir, 'mesh_cuboid', wall1_cuboid_ref, mesh_size, pp_tag)
547
548 os.system("mkdir -p ../out")
549
550 for s in p_mesh_fname:
551 print('\n\n')
552 print(s)
553 print("gmsh {}_{}.geo -3".format(s, pp_tag))
554 print('\n\n')
555 os.system("gmsh {}_{}.geo -3".format(s, pp_tag))
556 os.system("gmsh {}_{}.geo -3 -o {}_{}.vtk".format(s, pp_tag, s, pp_tag))
557
558
559
561inp_dir = './'
562pp_tag = 0
563if len(sys.argv) > 1:
564 pp_tag = int(sys.argv[1])
565
566create_input_file(inp_dir, pp_tag)
problem_setup.generate_cuboid_wall_gmsh_input
generate_cuboid_wall_gmsh_input(inp_dir, filename, cuboid, mesh_size, pp_tag)
Definition problem_setup.py:166
problem_setup.get_E
get_E(K, nu)
Definition problem_setup.py:82
problem_setup.get_eff_k
get_eff_k(k1, k2)
Definition problem_setup.py:89
problem_setup.generate_sphere_particle_gmsh_input
generate_sphere_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:784
problem_setup.particle_locations
particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R2, id_choices1, id_choices2, N1, N2, R_in, bar_rect, z_coord, add_orient=True)
Definition problem_setup.py:270
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:85