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])
74 sim_particles.append([1., R1, 2. * R1 + R2 + offset, 0., 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
83def particle_locations_orient(inp_dir, pp_tag, R1, R2, offset):
84 """Generate particle location data"""
85
86 sim_particles = []
87 sim_particles.append([0., R1, R1, 0., R1, 0.5*np.pi])
88 sim_particles.append([1., R1, 2. * R1 + R2 + offset, 0., R2, 0.])
89
90 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
91 inpf.write("i, x, y, z, r, o\n")
92 for p in sim_particles:
93 inpf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4], p[5]))
94
95 inpf.close()
96
97def rotate(p, theta, axis):
98
99 p_np = np.array(p)
100 axis_np = np.array(axis)
101
102 ct = np.cos(theta);
103 st = np.sin(theta);
104
105 # dot
106 p_dot_n = np.dot(p_np,axis_np)
107
108 # cross
109 n_cross_p = np.cross(axis_np, p_np)
110
111 return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p
112
113
114def get_ref_hex_points(center, radius, add_center = False):
115
116 # drum2d
117 #
118 # v3 v2
119 # + +
120 #
121 #
122 # + o +
123 # v4 x v1
124 #
125 # + +
126 # v5 v6
127 #
128 # Axis is a vector from x to v1
129 #
130
131 axis = [1., 0., 0.]
132
133 # center and radius
134 sim_Cx = center[0]
135 sim_Cy = center[1]
136 sim_Cz = center[2]
137 sim_radius = radius
138
139 rotate_axis = [0., 0., 1.]
140
141 points = []
142 if add_center:
143 points.append(center)
144
145 for i in range(6):
146 xi = rotate(axis, i*np.pi/3., rotate_axis)
147 points.append([center[i] + radius * xi[i] for i in range(3)])
148
149 return points
150
151
152def get_ref_drum_points(center, radius, width, add_center = False):
153
154 # drum2d
155 #
156 # v3 v2
157 # + +
158 #
159 #
160 # + o +
161 # v4 x v1
162 #
163 # + +
164 # v5 v6
165 #
166 # Axis is a vector from x to v1
167 #
168
169 axis = [1., 0., 0.]
170
171 # center and radius
172 sim_Cx = center[0]
173 sim_Cy = center[1]
174 sim_Cz = center[2]
175 sim_radius = radius
176
177 rotate_axis = [0., 0., 1.]
178
179 x1 = rotate(axis, np.pi/3., rotate_axis)
180 x2 = rotate(axis, -np.pi/3., rotate_axis)
181
182 points = []
183 if add_center:
184 points.append(center)
185
186 # v1
187 points.append([center[i] + width*0.5*axis[i] for i in range(3)])
188 # v2
189 points.append([center[i] + radius*x1[i] for i in range(3)])
190 # v3
191 points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)])
192 # v4
193 points.append([center[i] - width*0.5*axis[i] for i in range(3)])
194 # v5
195 v6 = [center[i] + radius*x2[i] for i in range(3)]
196 points.append([v6[i] - radius*axis[i] for i in range(3)])
197 # v6
198 points.append(v6)
199
200 return points
201
202
203def generate_hex_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
204
205 sim_inp_dir = str(inp_dir)
206
207 points = get_ref_hex_points(center, radius, True)
208
209 #
210 # create .geo file for gmsh
211 #
212 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
213 geof.write("cl__1 = 1;\n")
214 geof.write("Mesh.MshFileVersion = 2.2;\n")
215
216 #
217 # points
218 #
219 for i in range(7):
220 p = points[i]
221 sts = "Point({}) = ".format(i+1)
222 sts += "{"
223 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
224 sts += "};\n"
225 geof.write(sts);
226
227 #
228 # circlular arc
229 #
230 geof.write("Line(1) = {2, 3};\n")
231 geof.write("Line(2) = {3, 4};\n")
232 geof.write("Line(3) = {4, 5};\n")
233 geof.write("Line(4) = {5, 6};\n")
234 geof.write("Line(5) = {6, 7};\n")
235 geof.write("Line(6) = {7, 2};\n")
236
237 #
238 # surfaces
239 #
240 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
241
242 #
243 # plane surface
244 #
245 geof.write("Plane Surface(1) = {1};\n")
246
247 #
248 # physical surface
249 #
250 # tag = '"' + "a" + '"'
251 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
252
253 # # add center point to plane surface
254 geof.write("Point{1} In Surface {1};")
255
256 # close file
257 geof.close()
258
259
260def generate_drum_particle_gmsh_input(inp_dir, filename, center, radius, width, mesh_size, pp_tag):
261
262 sim_inp_dir = str(inp_dir)
263
264 # points
265 points = get_ref_drum_points(center, radius, width, True)
266
267 #
268 # create .geo file for gmsh
269 #
270 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
271 geof.write("cl__1 = 1;\n")
272 geof.write("Mesh.MshFileVersion = 2.2;\n")
273
274 #
275 # points
276 #
277 for i in range(7):
278 p = points[i]
279 sts = "Point({}) = ".format(i+1)
280 sts += "{"
281 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
282 sts += "};\n"
283 geof.write(sts);
284
285 #
286 # circlular arc
287 #
288 geof.write("Line(1) = {2, 3};\n")
289 geof.write("Line(2) = {3, 4};\n")
290 geof.write("Line(3) = {4, 5};\n")
291 geof.write("Line(4) = {5, 6};\n")
292 geof.write("Line(5) = {6, 7};\n")
293 geof.write("Line(6) = {7, 2};\n")
294
295 #
296 # surfaces
297 #
298 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
299
300 #
301 # plane surface
302 #
303 geof.write("Plane Surface(1) = {1};\n")
304
305 #
306 # physical surface
307 #
308 # tag = '"' + "a" + '"'
309 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
310
311 # # add center point to plane surface
312 geof.write("Point{1} In Surface {1};")
313
314 # close file
315 geof.close()
316
317
318
319def create_input_file(inp_dir, pp_tag):
320 """Generates input file for two-particle test"""
321
322 sim_inp_dir = str(inp_dir)
323
324
325 center = [0., 0., 0.]
326 R1 = 0.0015
327 R2 = 0.001
328 mesh_size = R1 / 10.
329 if R2 < R1:
330 mesh_size = R2 / 10.
331
332 horizon = 3. * mesh_size
333 particle_dist = 0.0004
334
335
336 final_time = 0.002
337 num_steps = 30000
338 # final_time = 0.00002
339 # num_steps = 2
340 num_outputs = 10
341 dt_out_n = num_steps / num_outputs
342 perform_out = True
343
344
345 poisson1 = 0.25
346 rho1 = 1200.
347 K1 = 2.16e+7
348 E1 = get_E(K1, poisson1)
349 G1 = get_G(E1, poisson1)
350 Gc1 = 50.
351
352 poisson2 = 0.25
353 rho2 = 1200.
354 K2 = 2.16e+7
355 E2 = get_E(K2, poisson2)
356 G2 = get_G(E2, poisson2)
357 Gc2 = 50.
358
359
362 R_contact_factor = 0.95
363
364 # Kn_V_max = 7.385158e+05
365 # Kn = np.power(Kn_V_max, 2)
366 # compute from bulk modulus
367
368 # from bulk modulus
369 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
370 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
371 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
372
373 beta_n_eps = 0.95
374 friction_coeff = 0.5
375 damping_active = False
376 friction_active = False
377 beta_n_factor = 100.
378
379
380 gravity_active = True
381 gravity = [0., -10., 0.]
382
383
384 free_fall_dist = particle_dist - horizon
385 free_fall_vel = [0., 0., 0.]
386 # free_fall_vel[1] = -np.sqrt(2. * np.abs(gravity[1]) * free_fall_dist)
387 free_fall_vel[1] = -1.
388
389
390 neigh_search_factor = 10.
391 neigh_search_interval = 40
392 neigh_search_criteria = "simple_all"
393
394
395
399 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
400 inpf.write("Model:\n")
401 inpf.write(" Dimension: 2\n")
402 inpf.write(" Discretization_Type:\n")
403 inpf.write(" Spatial: finite_difference\n")
404 inpf.write(" Time: central_difference\n")
405 inpf.write(" Final_Time: %4.6e\n" % (final_time))
406 inpf.write(" Time_Steps: %d\n" % (num_steps))
407
408 inpf.write("Policy:\n")
409 inpf.write(" Enable_PostProcessing: true\n")
410
411 #
412 # zone info
413 #
414 inpf.write("Zone:\n")
415 inpf.write(" Zones: 2\n")
416
417
418 inpf.write(" Zone_1:\n")
419 inpf.write(" Is_Wall: false\n")
420
421
422 inpf.write(" Zone_2:\n")
423 inpf.write(" Is_Wall: false\n")
424
425 #
426 # container info
427 #
428 inpf.write("Container:\n")
429 inpf.write(" Geometry:\n")
430 inpf.write(" Type: rectangle\n")
431 contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + particle_dist, 0.]
432 if R2 > R1:
433 contain_params[3] = 2.*R2
434 inpf.write(" Parameters: " + print_dbl_list(contain_params))
435
436 #
437 # particle info
438 #
439 inpf.write("Particle:\n")
440 inpf.write(" Test_Name: two_particle\n")
441 inpf.write(" Zone_1:\n")
442 inpf.write(" Type: drum2d\n")
443 drum_axis = [1., 0., 0.]
444 drum1_neck_width = 0.5*R1
445 p1_geom = [R1, drum1_neck_width, center[0], center[1], center[2], drum_axis[0], drum_axis[1], drum_axis[2]]
446 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
447 inpf.write(" Zone_2:\n")
448 inpf.write(" Type: hexagon\n")
449 hex_axis = [1., 0., 0.]
450 p2_geom = [R2, center[0], center[1], center[2], hex_axis[0], hex_axis[1], hex_axis[2]]
451 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
452
453 #
454 # particle generation
455 #
456 inpf.write("Particle_Generation:\n")
457 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
458 # specify that we also provide the orientation information in the file
459 inpf.write(" File_Data_Type: loc_rad_orient\n")
460
461 #
462 # Mesh info
463 #
464 inpf.write("Mesh:\n")
465
466
467 inpf.write(" Zone_1:\n")
468 inpf.write(" File: mesh_particle_1_" + str(pp_tag) + ".msh \n")
469
470
471 inpf.write(" Zone_2:\n")
472 inpf.write(" File: mesh_particle_2_" + str(pp_tag) + ".msh \n")
473
474 #
475 # Contact info
476 #
477 inpf.write("Contact:\n")
478
479
480 inpf.write(" Zone_11:\n")
481 # inpf.write(" Contact_Radius: %4.6e\n" % (R_contact))
482 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
483
484 if damping_active == False:
485 inpf.write(" Damping_On: false\n")
486 if friction_active == False:
487 inpf.write(" Friction_On: false\n")
488
489 inpf.write(" Kn: %4.6e\n" % (Kn_11))
490 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
491 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
492 inpf.write(" Kn_Factor: 1.0\n")
493 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
494
495
496 inpf.write(" Zone_12:\n")
497 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
498
499 if damping_active == False:
500 inpf.write(" Damping_On: false\n")
501 if friction_active == False:
502 inpf.write(" Friction_On: false\n")
503
504 inpf.write(" Kn: %4.6e\n" % (Kn_12))
505 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
506 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
507 inpf.write(" Kn_Factor: 1.0\n")
508 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
509
510
511 inpf.write(" Zone_22:\n")
512 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
513
514 if damping_active == False:
515 inpf.write(" Damping_On: false\n")
516 if friction_active == False:
517 inpf.write(" Friction_On: false\n")
518
519 inpf.write(" Kn: %4.6e\n" % (Kn_22))
520 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
521 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
522 inpf.write(" Kn_Factor: 1.0\n")
523 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
524
525 # Neighbor info
526 inpf.write("Neighbor:\n")
527 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
528 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
529 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
530
531 #
532 # Material info
533 #
534 inpf.write("Material:\n")
535
536
537 inpf.write(" Zone_1:\n")
538 inpf.write(" Type: PDState\n")
539 inpf.write(" Horizon: %4.6e\n" % (horizon))
540 inpf.write(" Density: %4.6e\n" % (rho1))
541 inpf.write(" Compute_From_Classical: true\n")
542 inpf.write(" K: %4.6e\n" % (K1))
543 inpf.write(" G: %4.6e\n" % (G1))
544 inpf.write(" Gc: %4.6e\n" % (Gc1))
545 inpf.write(" Influence_Function:\n")
546 inpf.write(" Type: 1\n")
547
548
549 inpf.write(" Zone_2:\n")
550 inpf.write(" Type: PDState\n")
551 inpf.write(" Horizon: %4.6e\n" % (horizon))
552 inpf.write(" Density: %4.6e\n" % (rho2))
553 inpf.write(" Compute_From_Classical: true\n")
554 inpf.write(" K: %4.6e\n" % (K2))
555 inpf.write(" G: %4.6e\n" % (G2))
556 inpf.write(" Gc: %4.6e\n" % (Gc2))
557 inpf.write(" Influence_Function:\n")
558 inpf.write(" Type: 1\n")
559
560 #
561 # Force
562 #
563 if gravity_active == True:
564 inpf.write("Force_BC:\n")
565 inpf.write(" Gravity: " + print_dbl_list(gravity))
566
567 #
568 # IC
569 #
570 inpf.write("IC:\n")
571 inpf.write(" Constant_Velocity:\n")
572 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
573 inpf.write(" Particle_List: [1]\n")
574
575 #
576 # Displacement
577 #
578 inpf.write("Displacement_BC:\n")
579 inpf.write(" Sets: 1\n")
580
581 inpf.write(" Set_1:\n")
582 inpf.write(" Particle_List: [0]\n")
583 inpf.write(" Direction: [1,2]\n")
584 inpf.write(" Time_Function:\n")
585 inpf.write(" Type: constant\n")
586 inpf.write(" Parameters:\n")
587 inpf.write(" - 0.0\n")
588 inpf.write(" Spatial_Function:\n")
589 inpf.write(" Type: constant\n")
590 inpf.write(" Zero_Displacement: true\n")
591
592 #
593 # Output info
594 #
595 inpf.write("Output:\n")
596 inpf.write(" Path: ../out/\n")
597 inpf.write(" Tags:\n")
598 inpf.write(" - Displacement\n")
599 inpf.write(" - Velocity\n")
600 inpf.write(" - Force\n")
601 inpf.write(" - Force_Density\n")
602 inpf.write(" - Damage_Z\n")
603 inpf.write(" - Damage\n")
604 inpf.write(" - Nodal_Volume\n")
605 inpf.write(" - Zone_ID\n")
606 inpf.write(" - Particle_ID\n")
607 inpf.write(" - Fixity\n")
608 inpf.write(" - Force_Fixity\n")
609 inpf.write(" - Contact_Data\n")
610 inpf.write(" - No_Fail_Node\n")
611 inpf.write(" - Boundary_Node_Flag\n")
612 inpf.write(" - Theta\n")
613 # inpf.write(" - Strain_Stress\n")
614 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
615 inpf.write(" Compress_Type: zlib\n")
616 inpf.write(" Perform_FE_Out: false\n")
617 if perform_out:
618 inpf.write(" Perform_Out: true\n")
619 else:
620 inpf.write(" Perform_Out: false\n")
621 inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))
622
623 inpf.write(" Debug: 1\n")
624 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
625 # inpf.write(" Output_Criteria: \n")
626 # inpf.write(" Type: max_particle_dist\n")
627 # inpf.write(" Parameters: [%4.6e]\n" % (2. * sim_h))
628
629 # close file
630 inpf.close()
631
632
633 # generate particle locations
634 # particle_locations(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
635 particle_locations_orient(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
636
637 p_mesh_fname = ['mesh_particle_1', 'mesh_particle_2']
638 # generate particle .geo file (large)
639 generate_drum_particle_gmsh_input(inp_dir, p_mesh_fname[0], center, R1, drum1_neck_width, mesh_size, pp_tag)
640 generate_hex_particle_gmsh_input(inp_dir, p_mesh_fname[1], center, R2, mesh_size, pp_tag)
641
642 os.system("mkdir -p ../out")
643
644 for s in p_mesh_fname:
645 print('\n\n')
646 print(s)
647 print("gmsh {}_{}.geo -2".format(s, pp_tag))
648 print('\n\n')
649 os.system("gmsh {}_{}.geo -2".format(s, pp_tag))
650 # os.system("gmsh {}_{}.geo -2 -o {}_{}.vtk".format(s, pp_tag, s, pp_tag))
651
652
653
655inp_dir = './'
656pp_tag = 0
657if len(sys.argv) > 1:
658 pp_tag = int(sys.argv[1])
659
660create_input_file(inp_dir, pp_tag)
problem_setup.generate_drum_particle_gmsh_input
generate_drum_particle_gmsh_input(inp_dir, filename, center, radius, width, mesh_size, pp_tag)
Definition problem_setup.py:260
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_hex_particle_gmsh_input
generate_hex_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag)
Definition problem_setup.py:535
problem_setup.create_input_file
create_input_file(inp_dir, pp_tag)
Definition problem_setup.py:793
problem_setup.particle_locations_orient
particle_locations_orient(inp_dir, pp_tag, R1, R2, offset)
Definition problem_setup.py:83
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
problem_setup.get_ref_drum_points
get_ref_drum_points(center, radius, width, add_center=False)
Definition problem_setup.py:236