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