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 # mw_dy1 = 0.032
577 mw_dy1 = dy1
578 mw_rect = [dx0-wpd, mw_dy1+wpd, dz0, dx1+wpd, mw_dy1+rwp, dz0]
579
580 # maximum distance to stop the simulation
581 max_dist_check = 20.*dx1
582
583
584 final_time = 0.06
585 num_steps = 600000
586 # final_time = 0.00002
587 # num_steps = 1
588 num_outputs = 400
589 dt_out_n = num_steps / num_outputs
590 test_dt_out_n = dt_out_n / 100
591 perform_out = True
592
593
594 poisson1 = 0.25
595 rho1 = 1200.
596 K1 = 2.16e+7
597 E1 = get_E(K1, poisson1)
598 G1 = get_G(E1, poisson1)
599 Gc1 = 50.
600
601 poisson2 = 0.25
602 rho2 = 1200.
603 K2 = 2.16e+7
604 E2 = get_E(K2, poisson2)
605 G2 = get_G(E2, poisson2)
606 Gc2 = 50.
607
608 # rigid wall
609 poisson3 = 0.25
610 rho3 = 1200.
611 K3 = 2.16e+7
612 E3 = get_E(K3, poisson3)
613 G3 = get_G(E3, poisson3)
614 Gc3 = 50.
615
616 # moving wall
617 poisson4 = 0.25
618 rho4 = 1200.
619 K4 = 2.16e+7
620 E4 = get_E(K4, poisson4)
621 G4 = get_G(E4, poisson4)
622 Gc4 = 50.
623
624
627 R_contact_factor = 0.95
628
629 # Kn_V_max = 7.385158e+05
630 # Kn = np.power(Kn_V_max, 2)
631 # compute from bulk modulus
632
633 # from bulk modulus
634 Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
635 Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
636 Kn_33 = 18. * get_eff_k(K3, K3) / (np.pi * np.power(horizon, 5))
637 Kn_44 = 18. * get_eff_k(K4, K4) / (np.pi * np.power(horizon, 5))
638 Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))
639 Kn_13 = 18. * get_eff_k(K1, K3) / (np.pi * np.power(horizon, 5))
640 Kn_14 = 18. * get_eff_k(K1, K4) / (np.pi * np.power(horizon, 5))
641 Kn_23 = 18. * get_eff_k(K2, K3) / (np.pi * np.power(horizon, 5))
642 Kn_24 = 18. * get_eff_k(K2, K4) / (np.pi * np.power(horizon, 5))
643 Kn_34 = 18. * get_eff_k(K3, K4) / (np.pi * np.power(horizon, 5))
644
645 beta_n_eps = 0.95
646 friction_coeff = 0.5
647 damping_active = True
648 friction_active = False
649 beta_n_factor = 100.
650
651
652 gravity_active = True
653 gravity = [0., -10., 0.]
654
655
656 mw_vy = -0.06
657
658
659
663 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
664 inpf.write("Model:\n")
665 inpf.write(" Dimension: 2\n")
666 inpf.write(" Discretization_Type:\n")
667 inpf.write(" Spatial: finite_difference\n")
668 inpf.write(" Time: central_difference\n")
669 inpf.write(" Final_Time: %4.6e\n" % (final_time))
670 inpf.write(" Time_Steps: %d\n" % (num_steps))
671
672 #
673 # container info
674 #
675 inpf.write("Container:\n")
676 inpf.write(" Geometry:\n")
677 inpf.write(" Type: rectangle\n")
678 inpf.write(" Parameters: " + print_dbl_list(rwo_rect))
679
680 #
681 # zone info
682 #
683 inpf.write("Zone:\n")
684 inpf.write(" Zones: 4\n")
685
686
687 inpf.write(" Zone_1:\n")
688 inpf.write(" Is_Wall: false\n")
689
690
691 inpf.write(" Zone_2:\n")
692 inpf.write(" Is_Wall: false\n")
693
694
695 inpf.write(" Zone_3:\n")
696 inpf.write(" Is_Wall: true\n")
697 inpf.write(" Type: rectangle\n")
698 inpf.write(" Parameters: " + print_dbl_list(rwo_rect))
699
700
701 inpf.write(" Zone_4:\n")
702 inpf.write(" Is_Wall: true\n")
703 inpf.write(" Type: rectangle\n")
704 inpf.write(" Parameters: " + print_dbl_list(mw_rect))
705
706 #
707 # particle info
708 #
709 inpf.write("Particle:\n")
710
711 inpf.write(" Test_Name: compressive_test\n")
712 inpf.write(" Compressive_Test:\n")
713 inpf.write(" Wall_Id: 1\n")
714 inpf.write(" Wall_Force_Direction: 2\n")
715
716 inpf.write(" Zone_1:\n")
717 inpf.write(" Type: circle\n")
718 p1_geom = [R1, center[0], center[1], center[2]]
719 inpf.write(" Parameters: " + print_dbl_list(p1_geom))
720 inpf.write(" Zone_2:\n")
721 inpf.write(" Type: hexagon\n")
722 hex_axis = [1., 0., 0.]
723 p2_geom = [R2, center[0], center[1], center[2], hex_axis[0], hex_axis[1], hex_axis[2]]
724 inpf.write(" Parameters: " + print_dbl_list(p2_geom))
725
726 #
727 # wall info
728 #
729 inpf.write("Wall:\n")
730 inpf.write(" Zone_3:\n")
731 inpf.write(" Type: flexible\n")
732 inpf.write(" All_Dofs_Constrained: true\n")
733 inpf.write(" Mesh: true\n")
734
735 inpf.write(" Zone_4:\n")
736 inpf.write(" Type: flexible\n")
737 inpf.write(" All_Dofs_Constrained: false\n")
738 inpf.write(" Mesh: true\n")
739
740 #
741 # particle generation
742 #
743 inpf.write("Particle_Generation:\n")
744 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
745 # specify that we also provide the orientation information in the file
746 inpf.write(" File_Data_Type: loc_rad_orient\n")
747
748 #
749 # Mesh info
750 #
751 inpf.write("Mesh:\n")
752
753
754 inpf.write(" Zone_1:\n")
755 inpf.write(" File: mesh_particle_1_" + str(pp_tag) + ".msh \n")
756
757
758 inpf.write(" Zone_2:\n")
759 inpf.write(" File: mesh_particle_2_" + str(pp_tag) + ".msh \n")
760
761
762 inpf.write(" Zone_3:\n")
763 inpf.write(" File: mesh_rigid_wall_" + str(pp_tag) + ".msh \n")
764
765
766 inpf.write(" Zone_4:\n")
767 inpf.write(" File: mesh_moving_wall_" + str(pp_tag) + ".msh \n")
768
769 #
770 # Contact info
771 #
772 inpf.write("Contact:\n")
773
774
775 inpf.write(" Zone_11:\n")
776 # inpf.write(" Contact_Radius: %4.6e\n" % (R_contact))
777 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
778
779 if damping_active == False:
780 inpf.write(" Damping_On: false\n")
781 if friction_active == False:
782 inpf.write(" Friction_On: false\n")
783
784 inpf.write(" Kn: %4.6e\n" % (Kn_11))
785 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
786 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
787 inpf.write(" Kn_Factor: 1.0\n")
788 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
789
790
791 inpf.write(" Zone_12:\n")
792 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
793
794 if damping_active == False:
795 inpf.write(" Damping_On: false\n")
796 if friction_active == False:
797 inpf.write(" Friction_On: false\n")
798
799 inpf.write(" Kn: %4.6e\n" % (Kn_12))
800 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
801 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
802 inpf.write(" Kn_Factor: 1.0\n")
803 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
804
805
806 inpf.write(" Zone_13:\n")
807 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
808
809 if damping_active == False:
810 inpf.write(" Damping_On: false\n")
811 if friction_active == False:
812 inpf.write(" Friction_On: false\n")
813
814 inpf.write(" Kn: %4.6e\n" % (Kn_13))
815 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
816 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
817 inpf.write(" Kn_Factor: 1.0\n")
818 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
819
820
821 inpf.write(" Zone_14:\n")
822 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
823
824 if damping_active == False:
825 inpf.write(" Damping_On: false\n")
826 if friction_active == False:
827 inpf.write(" Friction_On: false\n")
828
829 inpf.write(" Kn: %4.6e\n" % (Kn_14))
830 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
831 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
832 inpf.write(" Kn_Factor: 1.0\n")
833 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
834
835
836 inpf.write(" Zone_22:\n")
837 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
838
839 if damping_active == False:
840 inpf.write(" Damping_On: false\n")
841 if friction_active == False:
842 inpf.write(" Friction_On: false\n")
843
844 inpf.write(" Kn: %4.6e\n" % (Kn_22))
845 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
846 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
847 inpf.write(" Kn_Factor: 1.0\n")
848 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
849
850
851 inpf.write(" Zone_23:\n")
852 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
853
854 if damping_active == False:
855 inpf.write(" Damping_On: false\n")
856 if friction_active == False:
857 inpf.write(" Friction_On: false\n")
858
859 inpf.write(" Kn: %4.6e\n" % (Kn_23))
860 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
861 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
862 inpf.write(" Kn_Factor: 1.0\n")
863 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
864
865
866 inpf.write(" Zone_24:\n")
867 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
868
869 if damping_active == False:
870 inpf.write(" Damping_On: false\n")
871 if friction_active == False:
872 inpf.write(" Friction_On: false\n")
873
874 inpf.write(" Kn: %4.6e\n" % (Kn_24))
875 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
876 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
877 inpf.write(" Kn_Factor: 1.0\n")
878 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
879
880
881 inpf.write(" Zone_33:\n")
882 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
883
884 if damping_active == False:
885 inpf.write(" Damping_On: false\n")
886 if friction_active == False:
887 inpf.write(" Friction_On: false\n")
888
889 inpf.write(" Kn: %4.6e\n" % (Kn_33))
890 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
891 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
892 inpf.write(" Kn_Factor: 1.0\n")
893 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
894
895
896 inpf.write(" Zone_34:\n")
897 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
898
899 if damping_active == False:
900 inpf.write(" Damping_On: false\n")
901 if friction_active == False:
902 inpf.write(" Friction_On: false\n")
903
904 inpf.write(" Kn: %4.6e\n" % (Kn_34))
905 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
906 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
907 inpf.write(" Kn_Factor: 1.0\n")
908 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
909
910
911 inpf.write(" Zone_44:\n")
912 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
913
914 if damping_active == False:
915 inpf.write(" Damping_On: false\n")
916 if friction_active == False:
917 inpf.write(" Friction_On: false\n")
918
919 inpf.write(" Kn: %4.6e\n" % (Kn_44))
920 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
921 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
922 inpf.write(" Kn_Factor: 1.0\n")
923 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
924
925 #
926 # Neighbor info
927 #
928 inpf.write("Neighbor:\n")
929 inpf.write(" Update_Criteria: simple_all\n")
930 inpf.write(" Search_Factor: 5.0\n")
931
932 #
933 # Material info
934 #
935 inpf.write("Material:\n")
936
937
938 inpf.write(" Zone_1:\n")
939 inpf.write(" Type: PDState\n")
940 inpf.write(" Horizon: %4.6e\n" % (horizon))
941 inpf.write(" Density: %4.6e\n" % (rho1))
942 inpf.write(" Compute_From_Classical: true\n")
943 inpf.write(" K: %4.6e\n" % (K1))
944 inpf.write(" G: %4.6e\n" % (G1))
945 inpf.write(" Gc: %4.6e\n" % (Gc1))
946 inpf.write(" Influence_Function:\n")
947 inpf.write(" Type: 1\n")
948
949
950 inpf.write(" Zone_2:\n")
951 inpf.write(" Type: PDState\n")
952 inpf.write(" Horizon: %4.6e\n" % (horizon))
953 inpf.write(" Density: %4.6e\n" % (rho2))
954 inpf.write(" Compute_From_Classical: true\n")
955 inpf.write(" K: %4.6e\n" % (K2))
956 inpf.write(" G: %4.6e\n" % (G2))
957 inpf.write(" Gc: %4.6e\n" % (Gc2))
958 inpf.write(" Influence_Function:\n")
959 inpf.write(" Type: 1\n")
960
961
962 inpf.write(" Zone_3:\n")
963 inpf.write(" Type: PDState\n")
964 inpf.write(" Horizon: %4.6e\n" % (horizon))
965 inpf.write(" Density: %4.6e\n" % (rho3))
966 inpf.write(" Compute_From_Classical: true\n")
967 inpf.write(" K: %4.6e\n" % (K3))
968 inpf.write(" G: %4.6e\n" % (G3))
969 inpf.write(" Gc: %4.6e\n" % (Gc3))
970 inpf.write(" Influence_Function:\n")
971 inpf.write(" Type: 1\n")
972
973
974 inpf.write(" Zone_4:\n")
975 inpf.write(" Type: PDState\n")
976 inpf.write(" Horizon: %4.6e\n" % (horizon))
977 inpf.write(" Density: %4.6e\n" % (rho4))
978 inpf.write(" Compute_From_Classical: true\n")
979 inpf.write(" K: %4.6e\n" % (K4))
980 inpf.write(" G: %4.6e\n" % (G4))
981 inpf.write(" Gc: %4.6e\n" % (Gc4))
982 inpf.write(" Influence_Function:\n")
983 inpf.write(" Type: 1\n")
984
985 #
986 # Force
987 #
988 if gravity_active == True:
989 inpf.write("Force_BC:\n")
990 inpf.write(" Gravity: " + print_dbl_list(gravity))
991
992 #
993 # Displacement
994 #
995 inpf.write("Displacement_BC:\n")
996 inpf.write(" Sets: 2\n")
997
998 inpf.write(" Set_1:\n")
999 inpf.write(" Wall_List: [0]\n")
1000 inpf.write(" Direction: [1,2]\n")
1001 inpf.write(" Time_Function:\n")
1002 inpf.write(" Type: constant\n")
1003 inpf.write(" Parameters:\n")
1004 inpf.write(" - 0.0\n")
1005 inpf.write(" Spatial_Function:\n")
1006 inpf.write(" Type: constant\n")
1007 inpf.write(" Zero_Displacement: true\n")
1008
1009 inpf.write(" Set_2:\n")
1010 inpf.write(" Wall_List: [1]\n")
1011 inpf.write(" Direction: [2]\n")
1012 inpf.write(" Time_Function:\n")
1013 inpf.write(" Type: linear\n")
1014 inpf.write(" Parameters:\n")
1015 inpf.write(" - %4.6e\n" % (mw_vy))
1016 inpf.write(" Spatial_Function:\n")
1017 inpf.write(" Type: constant\n")
1018
1019 #
1020 # Output info
1021 #
1022 inpf.write("Output:\n")
1023 inpf.write(" Path: ../out/\n")
1024 inpf.write(" Tags:\n")
1025 inpf.write(" - Displacement\n")
1026 inpf.write(" - Velocity\n")
1027 inpf.write(" - Force\n")
1028 inpf.write(" - Force_Density\n")
1029 inpf.write(" - Damage_Z\n")
1030 inpf.write(" - Damage\n")
1031 inpf.write(" - Nodal_Volume\n")
1032 inpf.write(" - Zone_ID\n")
1033 inpf.write(" - Particle_ID\n")
1034 inpf.write(" - Fixity\n")
1035 inpf.write(" - Force_Fixity\n")
1036 inpf.write(" - Contact_Data\n")
1037 inpf.write(" - No_Fail_Node\n")
1038 inpf.write(" - Boundary_Node_Flag\n")
1039 # inpf.write(" - Particle_Locations\n")
1040 # inpf.write(" - Strain_Stress\n")
1041 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
1042 inpf.write(" Compress_Type: zlib\n")
1043 inpf.write(" Perform_FE_Out: false\n")
1044 if perform_out:
1045 inpf.write(" Perform_Out: true\n")
1046 else:
1047 inpf.write(" Perform_Out: false\n")
1048 inpf.write(" Test_Output_Interval: %d\n" % (test_dt_out_n))
1049
1050 inpf.write(" Debug: 1\n")
1051 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
1052 inpf.write(" Output_Criteria: \n")
1053 # inpf.write(" Type: max_particle_dist\n")
1054 # inpf.write(" Parameters: [%4.6e]\n" % (2. * sim_h))
1055 inpf.write(" Type: max_node_dist\n")
1056 inpf.write(" Parameters: [%4.6e]\n" % (max_dist_check))
1057
1058 #
1059 # restart
1060 #
1061 # inpf.write("Restart:\n")
1062 # inpf.write(" File: output_600\n")
1063 # inpf.write(" Step: 600000\n")
1064 # inpf.write(" Change_Reference_Free_Dofs: true\n")
1065 # inpf.write(" Tags:\n")
1066
1067 inpf.write("HPX:\n")
1068 inpf.write(" Partitions: 1\n")
1069
1070 # close file
1071 inpf.close()
1072
1073
1074 # generate particle locations
1075 re_generate_files = False
1076 if re_generate_files:
1077 particle_locations(inp_dir, pp_tag, R1, R2, [dx0, dy0, dz0, dx1, dy1, dz1], mesh_size, particle_padding)
1078
1079 # generate particle .geo file (large)
1080 generate_cir_particle_gmsh_input(inp_dir, 'mesh_particle_1', center, R1, mesh_size, pp_tag)
1081 generate_hex_particle_gmsh_input(inp_dir, 'mesh_particle_2', center, R2, mesh_size, pp_tag)
1082 generate_moving_wall_gmsh_input(inp_dir, 'mesh_moving_wall', mw_rect, mesh_size, pp_tag)
1083 generate_rigid_wall_gmsh_input(inp_dir, 'mesh_rigid_wall', rwo_rect, rwi_rect, mesh_size, pp_tag)
1084
1085
1087inp_dir = './'
1088pp_tag = 0
1089if len(sys.argv) > 1:
1090 pp_tag = int(sys.argv[1])
1091
1092create_input_file(inp_dir, pp_tag)