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 does_intersect_rect(p, r, particles, padding, rect):
59
60 # check intersection with rectangle
61 pr = [p[0] - r, p[1] - r, p[2], p[0] + r, p[1] + r, p[2]]
62
63 if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding:
64
65 # print('circle (xc = {}, r = {:5.3e}) intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
66
67 return True
68
69 # loop over particles
70 # for pi in particles:
71 # dx = [p[0] - pi[1], p[1] - pi[2], p[2] - pi[3]]
72 # rR = r + pi[4] + padding
73 # if np.linalg.norm(dx) < rR:
74 # return True
75
76 # print('circle (xc = {}, r = {:5.3e}) does not intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
77 return False
78
79
80def get_E(K, nu):
81 return 3. * K * (1. - 2. * nu)
82
83def get_G(E, nu):
84 return E / (2. * (1. + nu))
85
86
87def get_eff_k(k1, k2):
88 return 2. * k1 * k2 / (k1 + k2)
89
90
91def get_max(l):
92 l = np.array(l)
93 return l.max()
94
95
96def rotate(p, theta, axis):
97
98 p_np = np.array(p)
99 axis_np = np.array(axis)
100
101 ct = np.cos(theta);
102 st = np.sin(theta);
103
104 # dot
105 p_dot_n = np.dot(p_np,axis_np)
106
107 # cross
108 n_cross_p = np.cross(axis_np, p_np)
109
110 return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p
111
112
113def get_ref_hex_points(center, radius, add_center = False):
114
115 # drum2d
116 #
117 # v3 v2
118 # + +
119 #
120 #
121 # + o +
122 # v4 x v1
123 #
124 # + +
125 # v5 v6
126 #
127 # Axis is a vector from x to v1
128 #
129
130 axis = [1., 0., 0.]
131
132 # center and radius
133 sim_Cx = center[0]
134 sim_Cy = center[1]
135 sim_Cz = center[2]
136 sim_radius = radius
137
138 rotate_axis = [0., 0., 1.]
139
140 points = []
141 if add_center:
142 points.append(center)
143
144 for i in range(6):
145 xi = rotate(axis, i*np.pi/3., rotate_axis)
146 points.append([center[i] + radius * xi[i] for i in range(3)])
147
148 return points
149
150
151def particle_locations(inp_dir, pp_tag, R1, R2, rect, mesh_size, padding, add_orient = True):
152 """Generate particle location data"""
153
154 #
155 #
156 # Given rectangle, we want to fill it with particles of R1 and R2 radius
157 # in alternate fashion with padding between each particle
158 #
159 #
160 particles = []
161 points = []
162
163 # find how approximate number of rows and cols we can have
164 check_r = R1
165 if R1 > R2:
166 check_r = R2
167
168 rows = int((rect[3] - rect[0])/ check_r)
169 cols = int((rect[4] - rect[1])/ check_r)
170
171 rads = [R1, R2]
172
173 counter = 0
174 x_old = rect[0]
175 x_old_right = rect[3]
176 y_old = rect[1]
177 pad = padding
178
179 cx = 0.
180 cy = 0.
181 cz = 0.
182 r = 0.
183
184 row_rads = []
185 row_rads.append(R1)
186 row_rads.append(R2)
187
188 for i in range(rows):
189
190 if i > 0:
191 y_old = cy + get_max(row_rads)
192
193 # reset
194 row_rads = []
195 row_rads.append(R1)
196 row_rads.append(R2)
197
198 num_p_cols = 0
199 j = 0
200 while num_p_cols < cols - 1 and j < 500:
201
202 if j == 0:
203 # for odd row, take x_old as first and for even take as last
204 x_old = rect[0]
205 x_old_right = rect[3]
206
207 # ptype = counter % 2
208 ptype = 1
209 if np.random.uniform(0., 1.) < 0.7:
210 ptype = 0
211
212 r0 = rads[ptype]
213
214 # perturb radius
215 r = r0 + np.random.uniform(-0.1 * r0, 0.1 * r0)
216 # r = r0
217 row_rads.append(r)
218
219 # random horizontal perturbation by 10% of radius
220 rph = np.random.uniform(-0.1 * r0, 0.1 * r0)
221
222 cx0 = x_old + pad + r
223 cx = cx0 + rph
224 # for even row, we start from right edge of rectange
225 if i%2 == 0:
226 cx0 = x_old_right - pad - r
227 cx = cx0 - rph
228 cy = y_old + pad + r
229 cz = rect[2]
230
231 inters = does_intersect_rect([cx, cy, cz], r, particles, pad, rect)
232 inters_parts = does_intersect
233
234 if inters == False:
235 particles.append([float(ptype), cx, cy, cz, r])
236
237 # set x_old
238 if i % 2 == 0:
239 x_old_right = cx - r
240 else:
241 x_old = cx + r
242
243 counter += 1
244
245 num_p_cols += 1
246
247 j += 1
248
249
250 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
251 if add_orient:
252 inpf.write("i, x, y, z, r, o\n")
253 counter = 0
254 for p in particles:
255 inpf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4], np.random.uniform(0., 2.*np.pi)))
256 counter += 1
257 else:
258 inpf.write("i, x, y, z, r\n")
259 for p in particles:
260 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
261
262
263 inpf.close()
264
265 for p in particles:
266 points.append([p[1], p[2], p[3]])
267
268 points = np.array(points)
269 mesh = pv.PolyData(points)
270 pv.save_meshio('particle_locations_' + str(pp_tag) + '.vtu', mesh)
271
272 print('number of particles created = {}'.format(len(particles)))
273
274
275def generate_cir_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
276
277 sim_inp_dir = str(inp_dir)
278
279 # center and radius
280 sim_Cx = center[0]
281 sim_Cy = center[1]
282 sim_Cz = center[2]
283 sim_radius = radius
284
285 # mesh size
286 sim_h = mesh_size
287
288 #
289 # create .geo file for gmsh
290 #
291 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
292 geof.write("cl__1 = 1;\n")
293 geof.write("Mesh.MshFileVersion = 2.2;\n")
294
295 #
296 # points
297 #
298 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz, sim_h));
299 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx + sim_radius, sim_Cy, sim_Cz, sim_h))
300 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx - sim_radius, sim_Cy, sim_Cz, sim_h))
301 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy + sim_radius, sim_Cz, sim_h))
302 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy - sim_radius, sim_Cz, sim_h))
303
304 #
305 # circlular arc
306 #
307 geof.write("Circle(1) = {2, 1, 4};\n")
308 geof.write("Circle(2) = {4, 1, 3};\n")
309 geof.write("Circle(3) = {3, 1, 5};\n")
310 geof.write("Circle(4) = {5, 1, 2};\n")
311
312 #
313 # surfaces
314 #
315 geof.write("Line Loop(1) = {2, 3, 4, 1};\n")
316
317 #
318 # plane surface
319 #
320 geof.write("Plane Surface(1) = {1};\n")
321
322 #
323 # physical surface
324 #
325 # tag = '"' + "a" + '"'
326 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
327
328 # # add center point to plane surface
329 geof.write("Point{1} In Surface {1};")
330
331 # close file
332 geof.close()
333
334
335def generate_hex_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
336
337 sim_inp_dir = str(inp_dir)
338
339 points = get_ref_hex_points(center, radius, True)
340
341 #
342 # create .geo file for gmsh
343 #
344 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
345 geof.write("cl__1 = 1;\n")
346 geof.write("Mesh.MshFileVersion = 2.2;\n")
347
348 #
349 # points
350 #
351 for i in range(7):
352 p = points[i]
353 sts = "Point({}) = ".format(i+1)
354 sts += "{"
355 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
356 sts += "};\n"
357 geof.write(sts);
358
359 #
360 # circlular arc
361 #
362 geof.write("Line(1) = {2, 3};\n")
363 geof.write("Line(2) = {3, 4};\n")
364 geof.write("Line(3) = {4, 5};\n")
365 geof.write("Line(4) = {5, 6};\n")
366 geof.write("Line(5) = {6, 7};\n")
367 geof.write("Line(6) = {7, 2};\n")
368
369 #
370 # surfaces
371 #
372 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
373
374 #
375 # plane surface
376 #
377 geof.write("Plane Surface(1) = {1};\n")
378
379 #
380 # physical surface
381 #
382 # tag = '"' + "a" + '"'
383 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
384
385 # # add center point to plane surface
386 geof.write("Point{1} In Surface {1};")
387
388 # close file
389 geof.close()
390
391
392def generate_rigid_wall_gmsh_input(inp_dir, filename, outer_rect, inner_rect, mesh_size, pp_tag):
393
394 sim_inp_dir = str(inp_dir)
395
396 # outer rectangle
397 sim_Lx_out1 = outer_rect[0]
398 sim_Ly_out1 = outer_rect[1]
399 sim_Lz_out1 = outer_rect[2]
400 sim_Lx_out2 = outer_rect[3]
401 sim_Ly_out2 = outer_rect[4]
402 sim_Lz_out2 = outer_rect[5]
403
404 # inner rectangle
405 sim_Lx_in1 = inner_rect[0]
406 sim_Ly_in1 = inner_rect[1]
407 sim_Lz_in1 = inner_rect[2]
408 sim_Lx_in2 = inner_rect[3]
409 sim_Ly_in2 = inner_rect[4]
410 sim_Lz_in2 = inner_rect[5]
411
412 # mesh size
413 sim_h = mesh_size
414
415 #
416 # create .geo file for gmsh
417 #
418 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
419 geof.write("cl__1 = 1;\n")
420 geof.write("Mesh.MshFileVersion = 2.2;\n")
421
422 #
423 # L7 L3
424 # 4 + --- + 8 7 + --- + 3
425 # | | | |
426 # | | | |
427 # | | L6 L4 | |
428 # L8 | | | |
429 # | | | | L2
430 # | | L5 | |
431 # | + ------------------- + |
432 # | 5 6 |
433 # | |
434 # 1 + ------------------------------- + 2
435 # L1
436 #
437 #
438 # for point 8 and 7 --> choose y coordinate from outer rectangle and x
439 # coordinate from inner rectangle
440
441 #
442 # points
443 #
444 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h))
445 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
446 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
447 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
448
449 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_in1, sim_Lz_in1, sim_h))
450 geof.write("Point(6) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_in1, sim_Lz_in1, sim_h))
451 geof.write("Point(7) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_out2, sim_Lz_in1, sim_h))
452 geof.write("Point(8) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_out2, sim_Lz_in1, sim_h))
453
454 #
455 # lines
456 #
457 geof.write("Line(1) = {1, 2};\n")
458 geof.write("Line(2) = {2, 3};\n")
459 geof.write("Line(3) = {3, 7};\n")
460 geof.write("Line(4) = {7, 6};\n")
461 geof.write("Line(5) = {6, 5};\n")
462 geof.write("Line(6) = {5, 8};\n")
463 geof.write("Line(7) = {8, 4};\n")
464 geof.write("Line(8) = {4, 1};\n")
465
466 #
467 # surfaces
468 #
469 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8};\n")
470
471 #
472 # plane surface
473 #
474 geof.write("Plane Surface(1) = {1};\n")
475
476 #
477 # physical surface
478 #
479 tag = '"' + "a" + '"'
480 geof.write("Physical Surface(%s) = {1};\n" % (tag))
481
482 # close file
483 geof.close()
484
485
486def generate_moving_wall_gmsh_input(inp_dir, filename, rectangle, mesh_size, pp_tag):
487
488 sim_inp_dir = str(inp_dir)
489
490 # outer rectangle
491 sim_Lx_out1 = rectangle[0]
492 sim_Ly_out1 = rectangle[1]
493 sim_Lz_out1 = rectangle[2]
494 sim_Lx_out2 = rectangle[3]
495 sim_Ly_out2 = rectangle[4]
496 sim_Lz_out2 = rectangle[5]
497
498 # mesh size
499 sim_h = mesh_size
500
501 #
502 # create .geo file for gmsh
503 #
504 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
505 geof.write("cl__1 = 1;\n")
506 geof.write("Mesh.MshFileVersion = 2.2;\n")
507
508 #
509 # points
510 #
511 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h));
512 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
513 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
514 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
515
516 #
517 # lines
518 #
519 geof.write("Line(1) = {1, 2};\n")
520 geof.write("Line(2) = {2, 3};\n")
521 geof.write("Line(3) = {3, 4};\n")
522 geof.write("Line(4) = {4, 1};\n")
523
524 #
525 # surfaces
526 #
527 geof.write("Line Loop(1) = {1, 2, 3, 4};\n")
528
529 #
530 # plane surface
531 #
532 geof.write("Plane Surface(1) = {1};\n")
533
534 #
535 # physical surface
536 #
537 tag = '"' + "a" + '"'
538 geof.write("Physical Surface(%s) = {1};\n" % (tag))
539
540 # close file
541 geof.close()
542
543
544def create_input_file(inp_dir, pp_tag):
545 """Generates input file for two-particle test"""
546
547 sim_inp_dir = str(inp_dir)
548
549
550 dx0, dy0, dz0 = 0., 0., 0.
551 dx1, dy1, dz1 = 0.0685, 0.041, 0.
552
553
554 center = [0., 0., 0.]
555 R1 = 0.001
556 R2 = 0.001
557 mesh_size = R1 / 5.
558 if R2 < R1:
559 mesh_size = R2 / 5.
560
561 horizon = 3. * mesh_size
562 particle_dist = 0.001
563
564 #particle_padding = 2 * mesh_size
565 particle_padding = 0.3 * R1
566
567
568
569 # rigid wall
570 wpd = 1.1 * mesh_size # particle-wall initial distance
571 rwp = horizon + wpd # padding outside domain (thickness plust wpd)
572 rwo_rect = [dx0 - rwp, dy0 - rwp, dz0, dx1 + rwp, dy1 + rwp, dz0]
573 rwi_rect = [dx0 - wpd, dy0 - wpd, dz0, dx1 + wpd, dy1 + wpd, dz0]
574
575 # moving wall
576 # desired_y_begin = a - t*0.06 -> a = desired_y_begin + t*0.06
577 # mw_dy1 = dy1
578 mw_dy1 = 0.0346
579 mw_rect = [dx0-wpd, mw_dy1+wpd, dz0, dx1+wpd, mw_dy1+rwp, dz0]
580
581 # maximum distance to stop the simulation
582 max_dist_check = 20.*dx1
583
584
585 final_time = 0.15 # 0.06 from before + new time interval
586 num_steps = 1500000
587 # final_time = 0.00002
588 # num_steps = 1
589 num_outputs = 1000
590 dt_out_n = num_steps / num_outputs
591 test_dt_out_n = dt_out_n / 100
592 perform_out = True
593
594
595 poisson1 = 0.25
596 rho1 = 1200.
597 K1 = 2.16e+7
598 E1 = get_E(K1, poisson1)
599 G1 = get_G(E1, poisson1)
600 Gc1 = 50.
601
602 poisson2 = 0.25
603 rho2 = 1200.
604 K2 = 2.16e+7
605 E2 = get_E(K2, poisson2)
606 G2 = get_G(E2, poisson2)
607 Gc2 = 50.
608
609 # rigid wall
610 poisson3 = 0.25
611 rho3 = 1200.
612 K3 = 2.16e+7
613 E3 = get_E(K3, poisson3)
614 G3 = get_G(E3, poisson3)
615 Gc3 = 50.
616
617 # moving wall
618 poisson4 = 0.25
619 rho4 = 1200.
620 K4 = 2.16e+7
621 E4 = get_E(K4, poisson4)
622 G4 = get_G(E4, poisson4)
623 Gc4 = 50.
624
625
628 R_contact_factor = 0.95
629
630 # Kn_V_max = 7.385158e+05
631 # Kn = np.power(Kn_V_max, 2)
632 # compute from bulk modulus
633
634 # from bulk modulus
635 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
636 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
637 Kn_33 = 18. * get_eff_k(K3, K3) / (np.pi * np.power(horizon, 5))
638 Kn_44 = 18. * get_eff_k(K4, K4) / (np.pi * np.power(horizon, 5))
639 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
640 Kn_13 = 18. * get_eff_k(K1, K3) / (np.pi * np.power(horizon, 5))
641 Kn_14 = 18. * get_eff_k(K1, K4) / (np.pi * np.power(horizon, 5))
642 Kn_23 = 18. * get_eff_k(K2, K3) / (np.pi * np.power(horizon, 5))
643 Kn_24 = 18. * get_eff_k(K2, K4) / (np.pi * np.power(horizon, 5))
644 Kn_34 = 18. * get_eff_k(K3, K4) / (np.pi * np.power(horizon, 5))
645
646 beta_n_eps = 0.95
647 friction_coeff = 0.5
648 damping_active = True
649 friction_active = False
650 beta_n_factor = 100.
651
652
653 gravity_active = True
654 gravity = [0., -10., 0.]
655
656
657 mw_vy = -0.06
658
659
660
664 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
665 inpf.write("Model:\n")
666 inpf.write(" Dimension: 2\n")
667 inpf.write(" Discretization_Type:\n")
668 inpf.write(" Spatial: finite_difference\n")
669 inpf.write(" Time: central_difference\n")
670 inpf.write(" Final_Time: %4.6e\n" % (final_time))
671 inpf.write(" Time_Steps: %d\n" % (num_steps))
672
673 #
674 # container info
675 #
676 inpf.write("Container:\n")
677 inpf.write(" Geometry:\n")
678 inpf.write(" Type: rectangle\n")
679 inpf.write(" Parameters: " + print_dbl_list(rwo_rect))
680
681 #
682 # zone info
683 #
684 inpf.write("Zone:\n")
685 inpf.write(" Zones: 4\n")
686
687
688 inpf.write(" Zone_1:\n")
689 inpf.write(" Is_Wall: false\n")
690
691
692 inpf.write(" Zone_2:\n")
693 inpf.write(" Is_Wall: false\n")
694
695
696 inpf.write(" Zone_3:\n")
697 inpf.write(" Is_Wall: true\n")
698 inpf.write(" Type: rectangle\n")
699 inpf.write(" Parameters: " + print_dbl_list(rwo_rect))
700
701
702 inpf.write(" Zone_4:\n")
703 inpf.write(" Is_Wall: true\n")
704 inpf.write(" Type: rectangle\n")
705 inpf.write(" Parameters: " + print_dbl_list(mw_rect))
706
707 #
708 # particle info
709 #
710 inpf.write("Particle:\n")
711
712 inpf.write(" Test_Name: compressive_test\n")
713 inpf.write(" Compressive_Test:\n")
714 inpf.write(" Wall_Id: 1\n")
715 inpf.write(" Wall_Force_Direction: 2\n")
716
717 inpf.write(" Zone_1:\n")
718 inpf.write(" Type: circle\n")
719 p1_geom = [R1, center[0], center[1], center[2]]
720 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
721 inpf.write(" Zone_2:\n")
722 inpf.write(" Type: hexagon\n")
723 hex_axis = [1., 0., 0.]
724 p2_geom = [R2, center[0], center[1], center[2], hex_axis[0], hex_axis[1], hex_axis[2]]
725 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
726
727 #
728 # wall info
729 #
730 inpf.write("Wall:\n")
731 inpf.write(" Zone_3:\n")
732 inpf.write(" Type: flexible\n")
733 inpf.write(" All_Dofs_Constrained: true\n")
734 inpf.write(" Mesh: true\n")
735
736 inpf.write(" Zone_4:\n")
737 inpf.write(" Type: flexible\n")
738 inpf.write(" All_Dofs_Constrained: false\n")
739 inpf.write(" Mesh: true\n")
740
741 #
742 # particle generation
743 #
744 inpf.write("Particle_Generation:\n")
745 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
746 # specify that we also provide the orientation information in the file
747 inpf.write(" File_Data_Type: loc_rad_orient\n")
748
749 #
750 # Mesh info
751 #
752 inpf.write("Mesh:\n")
753
754
755 inpf.write(" Zone_1:\n")
756 inpf.write(" File: mesh_particle_1_" + str(pp_tag) + ".msh \n")
757
758
759 inpf.write(" Zone_2:\n")
760 inpf.write(" File: mesh_particle_2_" + str(pp_tag) + ".msh \n")
761
762
763 inpf.write(" Zone_3:\n")
764 inpf.write(" File: mesh_rigid_wall_" + str(pp_tag) + ".msh \n")
765
766
767 inpf.write(" Zone_4:\n")
768 inpf.write(" File: mesh_moving_wall_" + str(pp_tag) + ".msh \n")
769
770 #
771 # Contact info
772 #
773 inpf.write("Contact:\n")
774
775
776 inpf.write(" Zone_11:\n")
777 # inpf.write(" Contact_Radius: %4.6e\n" % (R_contact))
778 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
779
780 if damping_active == False:
781 inpf.write(" Damping_On: false\n")
782 if friction_active == False:
783 inpf.write(" Friction_On: false\n")
784
785 inpf.write(" Kn: %4.6e\n" % (Kn_11))
786 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
787 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
788 inpf.write(" Kn_Factor: 1.0\n")
789 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
790
791
792 inpf.write(" Zone_12:\n")
793 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
794
795 if damping_active == False:
796 inpf.write(" Damping_On: false\n")
797 if friction_active == False:
798 inpf.write(" Friction_On: false\n")
799
800 inpf.write(" Kn: %4.6e\n" % (Kn_12))
801 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
802 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
803 inpf.write(" Kn_Factor: 1.0\n")
804 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
805
806
807 inpf.write(" Zone_13:\n")
808 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
809
810 if damping_active == False:
811 inpf.write(" Damping_On: false\n")
812 if friction_active == False:
813 inpf.write(" Friction_On: false\n")
814
815 inpf.write(" Kn: %4.6e\n" % (Kn_13))
816 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
817 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
818 inpf.write(" Kn_Factor: 1.0\n")
819 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
820
821
822 inpf.write(" Zone_14:\n")
823 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
824
825 if damping_active == False:
826 inpf.write(" Damping_On: false\n")
827 if friction_active == False:
828 inpf.write(" Friction_On: false\n")
829
830 inpf.write(" Kn: %4.6e\n" % (Kn_14))
831 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
832 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
833 inpf.write(" Kn_Factor: 1.0\n")
834 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
835
836
837 inpf.write(" Zone_22:\n")
838 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
839
840 if damping_active == False:
841 inpf.write(" Damping_On: false\n")
842 if friction_active == False:
843 inpf.write(" Friction_On: false\n")
844
845 inpf.write(" Kn: %4.6e\n" % (Kn_22))
846 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
847 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
848 inpf.write(" Kn_Factor: 1.0\n")
849 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
850
851
852 inpf.write(" Zone_23:\n")
853 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
854
855 if damping_active == False:
856 inpf.write(" Damping_On: false\n")
857 if friction_active == False:
858 inpf.write(" Friction_On: false\n")
859
860 inpf.write(" Kn: %4.6e\n" % (Kn_23))
861 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
862 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
863 inpf.write(" Kn_Factor: 1.0\n")
864 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
865
866
867 inpf.write(" Zone_24:\n")
868 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
869
870 if damping_active == False:
871 inpf.write(" Damping_On: false\n")
872 if friction_active == False:
873 inpf.write(" Friction_On: false\n")
874
875 inpf.write(" Kn: %4.6e\n" % (Kn_24))
876 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
877 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
878 inpf.write(" Kn_Factor: 1.0\n")
879 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
880
881
882 inpf.write(" Zone_33:\n")
883 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
884
885 if damping_active == False:
886 inpf.write(" Damping_On: false\n")
887 if friction_active == False:
888 inpf.write(" Friction_On: false\n")
889
890 inpf.write(" Kn: %4.6e\n" % (Kn_33))
891 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
892 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
893 inpf.write(" Kn_Factor: 1.0\n")
894 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
895
896
897 inpf.write(" Zone_34:\n")
898 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
899
900 if damping_active == False:
901 inpf.write(" Damping_On: false\n")
902 if friction_active == False:
903 inpf.write(" Friction_On: false\n")
904
905 inpf.write(" Kn: %4.6e\n" % (Kn_34))
906 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
907 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
908 inpf.write(" Kn_Factor: 1.0\n")
909 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
910
911
912 inpf.write(" Zone_44:\n")
913 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
914
915 if damping_active == False:
916 inpf.write(" Damping_On: false\n")
917 if friction_active == False:
918 inpf.write(" Friction_On: false\n")
919
920 inpf.write(" Kn: %4.6e\n" % (Kn_44))
921 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
922 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
923 inpf.write(" Kn_Factor: 1.0\n")
924 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
925
926 #
927 # Neighbor info
928 #
929 inpf.write("Neighbor:\n")
930 inpf.write(" Update_Criteria: simple_all\n")
931 inpf.write(" Search_Factor: 5.0\n")
932
933 #
934 # Material info
935 #
936 inpf.write("Material:\n")
937
938
939 inpf.write(" Zone_1:\n")
940 inpf.write(" Type: PDState\n")
941 inpf.write(" Horizon: %4.6e\n" % (horizon))
942 inpf.write(" Density: %4.6e\n" % (rho1))
943 inpf.write(" Compute_From_Classical: true\n")
944 inpf.write(" K: %4.6e\n" % (K1))
945 inpf.write(" G: %4.6e\n" % (G1))
946 inpf.write(" Gc: %4.6e\n" % (Gc1))
947 inpf.write(" Influence_Function:\n")
948 inpf.write(" Type: 1\n")
949
950
951 inpf.write(" Zone_2:\n")
952 inpf.write(" Type: PDState\n")
953 inpf.write(" Horizon: %4.6e\n" % (horizon))
954 inpf.write(" Density: %4.6e\n" % (rho2))
955 inpf.write(" Compute_From_Classical: true\n")
956 inpf.write(" K: %4.6e\n" % (K2))
957 inpf.write(" G: %4.6e\n" % (G2))
958 inpf.write(" Gc: %4.6e\n" % (Gc2))
959 inpf.write(" Influence_Function:\n")
960 inpf.write(" Type: 1\n")
961
962
963 inpf.write(" Zone_3:\n")
964 inpf.write(" Type: PDState\n")
965 inpf.write(" Horizon: %4.6e\n" % (horizon))
966 inpf.write(" Density: %4.6e\n" % (rho3))
967 inpf.write(" Compute_From_Classical: true\n")
968 inpf.write(" K: %4.6e\n" % (K3))
969 inpf.write(" G: %4.6e\n" % (G3))
970 inpf.write(" Gc: %4.6e\n" % (Gc3))
971 inpf.write(" Influence_Function:\n")
972 inpf.write(" Type: 1\n")
973
974
975 inpf.write(" Zone_4:\n")
976 inpf.write(" Type: PDState\n")
977 inpf.write(" Horizon: %4.6e\n" % (horizon))
978 inpf.write(" Density: %4.6e\n" % (rho4))
979 inpf.write(" Compute_From_Classical: true\n")
980 inpf.write(" K: %4.6e\n" % (K4))
981 inpf.write(" G: %4.6e\n" % (G4))
982 inpf.write(" Gc: %4.6e\n" % (Gc4))
983 inpf.write(" Influence_Function:\n")
984 inpf.write(" Type: 1\n")
985
986 #
987 # Force
988 #
989 if gravity_active == True:
990 inpf.write("Force_BC:\n")
991 inpf.write(" Gravity: " + print_dbl_list(gravity))
992
993 #
994 # Displacement
995 #
996 inpf.write("Displacement_BC:\n")
997 inpf.write(" Sets: 2\n")
998
999 inpf.write(" Set_1:\n")
1000 inpf.write(" Wall_List: [0]\n")
1001 inpf.write(" Direction: [1,2]\n")
1002 inpf.write(" Time_Function:\n")
1003 inpf.write(" Type: constant\n")
1004 inpf.write(" Parameters:\n")
1005 inpf.write(" - 0.0\n")
1006 inpf.write(" Spatial_Function:\n")
1007 inpf.write(" Type: constant\n")
1008 inpf.write(" Zero_Displacement: true\n")
1009
1010 inpf.write(" Set_2:\n")
1011 inpf.write(" Wall_List: [1]\n")
1012 inpf.write(" Direction: [2]\n")
1013 inpf.write(" Time_Function:\n")
1014 inpf.write(" Type: linear\n")
1015 inpf.write(" Parameters:\n")
1016 inpf.write(" - %4.6e\n" % (mw_vy))
1017 inpf.write(" Spatial_Function:\n")
1018 inpf.write(" Type: constant\n")
1019
1020 #
1021 # Output info
1022 #
1023 inpf.write("Output:\n")
1024 inpf.write(" Path: ../out/\n")
1025 inpf.write(" Tags:\n")
1026 inpf.write(" - Displacement\n")
1027 inpf.write(" - Velocity\n")
1028 inpf.write(" - Force\n")
1029 inpf.write(" - Force_Density\n")
1030 inpf.write(" - Damage_Z\n")
1031 inpf.write(" - Damage\n")
1032 inpf.write(" - Nodal_Volume\n")
1033 inpf.write(" - Zone_ID\n")
1034 inpf.write(" - Particle_ID\n")
1035 inpf.write(" - Fixity\n")
1036 inpf.write(" - Force_Fixity\n")
1037 # inpf.write(" - Contact_Data\n")
1038 inpf.write(" - No_Fail_Node\n")
1039 inpf.write(" - Boundary_Node_Flag\n")
1040 # inpf.write(" - Particle_Locations\n")
1041 # inpf.write(" - Strain_Stress\n")
1042 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
1043 inpf.write(" Compress_Type: zlib\n")
1044 inpf.write(" Perform_FE_Out: false\n")
1045 if perform_out:
1046 inpf.write(" Perform_Out: true\n")
1047 else:
1048 inpf.write(" Perform_Out: false\n")
1049 inpf.write(" Test_Output_Interval: %d\n" % (test_dt_out_n))
1050
1051 inpf.write(" Debug: 1\n")
1052 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
1053 inpf.write(" Output_Criteria: \n")
1054 # inpf.write(" Type: max_particle_dist\n")
1055 # inpf.write(" Parameters: [%4.6e]\n" % (2. * sim_h))
1056 inpf.write(" Type: max_node_dist\n")
1057 inpf.write(" Parameters: [%4.6e]\n" % (max_dist_check))
1058
1059 #
1060 # restart
1061 #
1062 inpf.write("Restart:\n")
1063 inpf.write(" File: output_400\n")
1064 inpf.write(" Step: 600000\n")
1065
1066 inpf.write("HPX:\n")
1067 inpf.write(" Partitions: 1\n")
1068
1069 # close file
1070 inpf.close()
1071
1072
1073 # generate particle locations
1074 re_generate_files = False
1075 if re_generate_files:
1076 particle_locations(inp_dir, pp_tag, R1, R2, [dx0, dy0, dz0, dx1, dy1, dz1], mesh_size, particle_padding)
1077
1078 # generate particle .geo file (large)
1079 generate_cir_particle_gmsh_input(inp_dir, 'mesh_particle_1', center, R1, mesh_size, pp_tag)
1080 generate_hex_particle_gmsh_input(inp_dir, 'mesh_particle_2', center, R2, mesh_size, pp_tag)
1081 generate_moving_wall_gmsh_input(inp_dir, 'mesh_moving_wall', mw_rect, mesh_size, pp_tag)
1082 generate_rigid_wall_gmsh_input(inp_dir, 'mesh_rigid_wall', rwo_rect, rwi_rect, mesh_size, pp_tag)
1083
1084
1086inp_dir = './'
1087pp_tag = 0
1088if len(sys.argv) > 1:
1089 pp_tag = int(sys.argv[1])
1090
1091create_input_file(inp_dir, pp_tag)