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+1.*R2*np.sin(np.pi/3.), 2. * R1 + R2 + offset, 0., R2, 0.5*np.pi])
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_drum_points(center, radius, width, 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 x1 = rotate(axis, np.pi/3., rotate_axis)
142 x2 = rotate(axis, -np.pi/3., rotate_axis)
143
144 points = []
145 if add_center:
146 points.append(center)
147
148 # v1
149 points.append([center[i] + width*0.5*axis[i] for i in range(3)])
150 # v2
151 points.append([center[i] + radius*x1[i] for i in range(3)])
152 # v3
153 points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)])
154 # v4
155 points.append([center[i] - width*0.5*axis[i] for i in range(3)])
156 # v5
157 v6 = [center[i] + radius*x2[i] for i in range(3)]
158 points.append([v6[i] - radius*axis[i] for i in range(3)])
159 # v6
160 points.append(v6)
161
162 return points
163
164
165def generate_particle_gmsh_input(inp_dir, filename, center, radius, width, mesh_size, pp_tag):
166
167 sim_inp_dir = str(inp_dir)
168
169 # points
170 points = get_ref_drum_points(center, radius, width, True)
171
172 #
173 # create .geo file for gmsh
174 #
175 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
176 geof.write("cl__1 = 1;\n")
177 geof.write("Mesh.MshFileVersion = 2.2;\n")
178
179 #
180 # points
181 #
182 for i in range(7):
183 p = points[i]
184 sts = "Point({}) = ".format(i+1)
185 sts += "{"
186 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
187 sts += "};\n"
188 geof.write(sts);
189
190 #
191 # circlular arc
192 #
193 geof.write("Line(1) = {2, 3};\n")
194 geof.write("Line(2) = {3, 4};\n")
195 geof.write("Line(3) = {4, 5};\n")
196 geof.write("Line(4) = {5, 6};\n")
197 geof.write("Line(5) = {6, 7};\n")
198 geof.write("Line(6) = {7, 2};\n")
199
200 #
201 # surfaces
202 #
203 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
204
205 #
206 # plane surface
207 #
208 geof.write("Plane Surface(1) = {1};\n")
209
210 #
211 # physical surface
212 #
213 # tag = '"' + "a" + '"'
214 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
215
216 # # add center point to plane surface
217 geof.write("Point{1} In Surface {1};")
218
219 # close file
220 geof.close()
221
222
223
224def create_input_file(inp_dir, pp_tag):
225 """Generates input file for two-particle test"""
226
227 sim_inp_dir = str(inp_dir)
228
229
230 center = [0., 0., 0.]
231 R1 = 0.0015
232 R2 = 0.001
233 mesh_size = R1 / 10.
234 if R2 < R1:
235 mesh_size = R2 / 10.
236
237 horizon = 3. * mesh_size
238 particle_dist = 0.001
239
240
241 final_time = 0.01
242 num_steps = 50000
243 # final_time = 0.00002
244 # num_steps = 2
245 num_outputs = 10
246 dt_out_n = num_steps / num_outputs
247 perform_out = True
248
249
250 poisson1 = 0.25
251 rho1 = 1200.
252 K1 = 2.16e+7
253 E1 = get_E(K1, poisson1)
254 G1 = get_G(E1, poisson1)
255 Gc1 = 50.
256
257 poisson2 = 0.25
258 rho2 = 1200.
259 K2 = 2.16e+7
260 E2 = get_E(K2, poisson2)
261 G2 = get_G(E2, poisson2)
262 Gc2 = 50.
263
264
267 R_contact_factor = 0.95
268
269 # Kn_V_max = 7.385158e+05
270 # Kn = np.power(Kn_V_max, 2)
271 # compute from bulk modulus
272
273 # from bulk modulus
274 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
275 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
276 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
277
278 beta_n_eps = 0.95
279 friction_coeff = 0.5
280 damping_active = True
281 friction_active = False
282 beta_n_factor = 100.
283
284
285 gravity_active = True
286 gravity = [0., -10., 0.]
287
288
289 free_fall_dist = particle_dist - horizon
290 free_fall_vel = [0., 0., 0.]
291 # free_fall_vel[1] = -np.sqrt(2. * np.abs(gravity[1]) * free_fall_dist)
292 free_fall_vel[1] = -8.
293
294
295 neigh_search_factor = 10.
296 neigh_search_interval = 40
297 neigh_search_criteria = "simple_all"
298
299
300
304 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
305 inpf.write("Model:\n")
306 inpf.write(" Dimension: 2\n")
307 inpf.write(" Discretization_Type:\n")
308 inpf.write(" Spatial: finite_difference\n")
309 inpf.write(" Time: central_difference\n")
310 inpf.write(" Final_Time: %4.6e\n" % (final_time))
311 inpf.write(" Time_Steps: %d\n" % (num_steps))
312
313 inpf.write("Policy:\n")
314 inpf.write(" Enable_PostProcessing: true\n")
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 # container info
332 #
333 inpf.write("Container:\n")
334 inpf.write(" Geometry:\n")
335 inpf.write(" Type: rectangle\n")
336 contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + particle_dist, 0.]
337 if R2 > R1:
338 contain_params[3] = 2.*R2
339 inpf.write(" Parameters: " + print_dbl_list(contain_params))
340
341 #
342 # particle info
343 #
344 inpf.write("Particle:\n")
345 inpf.write(" Test_Name: two_particle\n")
346 inpf.write(" Zone_1:\n")
347 inpf.write(" Type: drum2d\n")
348 drum_axis = [1., 0., 0.]
349 drum1_neck_width = 0.5*R1
350 p1_geom = [R1, drum1_neck_width, center[0], center[1], center[2], drum_axis[0], drum_axis[1], drum_axis[2]]
351 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
352 inpf.write(" Zone_2:\n")
353 inpf.write(" Type: drum2d\n")
354
355 drum2_neck_width = 0.5*R2
356 p2_geom = [R2, drum2_neck_width, center[0], center[1], center[2], drum_axis[0], drum_axis[1], drum_axis[2]]
357 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
358
359 #
360 # particle generation
361 #
362 inpf.write("Particle_Generation:\n")
363 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
364 # specify that we also provide the orientation information in the file
365 inpf.write(" File_Data_Type: loc_rad_orient\n")
366
367 #
368 # Mesh info
369 #
370 inpf.write("Mesh:\n")
371
372
373 inpf.write(" Zone_1:\n")
374 inpf.write(" File: mesh_drum_1_" + str(pp_tag) + ".msh \n")
375
376
377 inpf.write(" Zone_2:\n")
378 inpf.write(" File: mesh_drum_2_" + str(pp_tag) + ".msh \n")
379
380 #
381 # Contact info
382 #
383 inpf.write("Contact:\n")
384
385
386 inpf.write(" Zone_11:\n")
387 # inpf.write(" Contact_Radius: %4.6e\n" % (R_contact))
388 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
389
390 if damping_active == False:
391 inpf.write(" Damping_On: false\n")
392 if friction_active == False:
393 inpf.write(" Friction_On: false\n")
394
395 inpf.write(" Kn: %4.6e\n" % (Kn_11))
396 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
397 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
398 inpf.write(" Kn_Factor: 1.0\n")
399 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
400
401
402 inpf.write(" Zone_12:\n")
403 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
404
405 if damping_active == False:
406 inpf.write(" Damping_On: false\n")
407 if friction_active == False:
408 inpf.write(" Friction_On: false\n")
409
410 inpf.write(" Kn: %4.6e\n" % (Kn_12))
411 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
412 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
413 inpf.write(" Kn_Factor: 1.0\n")
414 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
415
416
417 inpf.write(" Zone_22:\n")
418 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
419
420 if damping_active == False:
421 inpf.write(" Damping_On: false\n")
422 if friction_active == False:
423 inpf.write(" Friction_On: false\n")
424
425 inpf.write(" Kn: %4.6e\n" % (Kn_22))
426 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
427 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
428 inpf.write(" Kn_Factor: 1.0\n")
429 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
430
431 # Neighbor info
432 inpf.write("Neighbor:\n")
433 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
434 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
435 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
436
437 #
438 # Material info
439 #
440 inpf.write("Material:\n")
441
442
443 inpf.write(" Zone_1:\n")
444 inpf.write(" Type: PDState\n")
445 inpf.write(" Horizon: %4.6e\n" % (horizon))
446 inpf.write(" Density: %4.6e\n" % (rho1))
447 inpf.write(" Compute_From_Classical: true\n")
448 inpf.write(" K: %4.6e\n" % (K1))
449 inpf.write(" G: %4.6e\n" % (G1))
450 inpf.write(" Gc: %4.6e\n" % (Gc1))
451 inpf.write(" Influence_Function:\n")
452 inpf.write(" Type: 1\n")
453
454
455 inpf.write(" Zone_2:\n")
456 inpf.write(" Type: PDState\n")
457 inpf.write(" Horizon: %4.6e\n" % (horizon))
458 inpf.write(" Density: %4.6e\n" % (rho2))
459 inpf.write(" Compute_From_Classical: true\n")
460 inpf.write(" K: %4.6e\n" % (K2))
461 inpf.write(" G: %4.6e\n" % (G2))
462 inpf.write(" Gc: %4.6e\n" % (Gc2))
463 inpf.write(" Influence_Function:\n")
464 inpf.write(" Type: 1\n")
465
466 #
467 # Force
468 #
469 if gravity_active == True:
470 inpf.write("Force_BC:\n")
471 inpf.write(" Gravity: " + print_dbl_list(gravity))
472
473 #
474 # IC
475 #
476 inpf.write("IC:\n")
477 inpf.write(" Constant_Velocity:\n")
478 inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
479 inpf.write(" Particle_List: [1]\n")
480
481 #
482 # Displacement
483 #
484 inpf.write("Displacement_BC:\n")
485 inpf.write(" Sets: 1\n")
486
487 inpf.write(" Set_1:\n")
488 inpf.write(" Particle_List: [0]\n")
489 inpf.write(" Direction: [1,2]\n")
490 inpf.write(" Time_Function:\n")
491 inpf.write(" Type: constant\n")
492 inpf.write(" Parameters:\n")
493 inpf.write(" - 0.0\n")
494 inpf.write(" Spatial_Function:\n")
495 inpf.write(" Type: constant\n")
496 inpf.write(" Zero_Displacement: true\n")
497
498 #
499 # Output info
500 #
501 inpf.write("Output:\n")
502 inpf.write(" Path: ../out/\n")
503 inpf.write(" Tags:\n")
504 inpf.write(" - Displacement\n")
505 inpf.write(" - Velocity\n")
506 inpf.write(" - Force\n")
507 inpf.write(" - Force_Density\n")
508 inpf.write(" - Damage_Z\n")
509 inpf.write(" - Damage\n")
510 inpf.write(" - Nodal_Volume\n")
511 inpf.write(" - Zone_ID\n")
512 inpf.write(" - Particle_ID\n")
513 inpf.write(" - Fixity\n")
514 inpf.write(" - Force_Fixity\n")
515 inpf.write(" - Contact_Nodes\n")
516 inpf.write(" - No_Fail_Node\n")
517 inpf.write(" - Boundary_Node_Flag\n")
518 inpf.write(" - Theta\n")
519 # inpf.write(" - Strain_Stress\n")
520 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
521 inpf.write(" Compress_Type: zlib\n")
522 inpf.write(" Perform_FE_Out: false\n")
523 if perform_out:
524 inpf.write(" Perform_Out: true\n")
525 else:
526 inpf.write(" Perform_Out: false\n")
527 inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))
528
529 inpf.write(" Debug: 1\n")
530 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
531 # inpf.write(" Output_Criteria: \n")
532 # inpf.write(" Type: max_particle_dist\n")
533 # inpf.write(" Parameters: [%4.6e]\n" % (2. * sim_h))
534
535 # close file
536 inpf.close()
537
538
539 # generate particle locations
540 # particle_locations(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
541 particle_locations_orient(inp_dir, pp_tag, R1, R2, particle_dist - free_fall_dist)
542
543 p_mesh_fname = ['mesh_drum_1', 'mesh_drum_2']
544 # generate particle .geo file (large)
545 generate_particle_gmsh_input(inp_dir, p_mesh_fname[0], center, R1, drum1_neck_width, mesh_size, pp_tag)
546 generate_particle_gmsh_input(inp_dir, p_mesh_fname[1], center, R2, drum2_neck_width, 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 -2".format(s, pp_tag))
554 print('\n\n')
555 os.system("gmsh {}_{}.geo -2".format(s, pp_tag))
556 # os.system("gmsh {}_{}.geo -2 -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)