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
3# import csv
4import sys
5# import pyvista as pv
6
7def print_bool(arg, prefix = ""):
8
9 str = prefix
10 if arg == True:
11 str += "True\n"
12 else:
13 str += "False\n"
14 return str
15
16def print_dbl(arg, prefix = ""):
17
18 str = prefix + "%4.6e\n" % (arg)
19 return str
20
21def print_int(arg, prefix = ""):
22 str = prefix + "%d\n" % (arg)
23 return str
24
25def print_dbl_list(arg, prefix = ""):
26 str = prefix + "["
27 N = len(arg)
28 for i in range(N):
29 str += "%4.6e" % (arg[i])
30 if i < N - 1:
31 str += ", "
32 else:
33 str += "]\n"
34
35 return str
36
37def print_int_list(arg, prefix = ""):
38 str = prefix + "["
39 N = len(arg)
40 for i in range(N):
41 str += "%d" % (arg[i])
42 if i < N - 1:
43 str += ", "
44 else:
45 str += "]\n"
46
47 return str
48
49def serialize_matrix_list(p):
50 s = []
51 for q in p:
52 for w in q:
53 s.append(w)
54
55 return s
56
57
58def write_point_geo(geof, p_id, x, h):
59 geof.write("Point(%d) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (p_id, x[0], x[1], x[2], h))
60 return p_id + 1
61
62def write_line_geo(geof, l_id, p1_id, p2_id):
63 geof.write("Line(%d) = {%d, %d};\n" % (l_id, p1_id, p2_id))
64 return l_id + 1
65
66def write_cir_line_geo(geof, l_id, p1_id, p2_id, p3_id):
67 geof.write("Circle(%d) = {%d, %d, %d};\n" % (l_id, p1_id, p2_id, p3_id))
68 return l_id + 1
69
70def write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, zone_string, Kn):
71 inpf.write(" Zone_%s:\n" % (zone_string))
72 inpf.write(" Contact_Radius_Factor: %4.6e\n" % (R_contact_factor))
73
74 if damping_active == False:
75 inpf.write(" Damping_On: false\n")
76 if friction_active == False:
77 inpf.write(" Friction_On: false\n")
78
79 inpf.write(" Kn: %4.6e\n" % (Kn))
80 inpf.write(" Epsilon: %4.6e\n" % (beta_n_eps))
81 inpf.write(" Friction_Coeff: %4.6e\n" % (friction_coeff))
82 inpf.write(" Kn_Factor: %4.6e\n" % (Kn_factor))
83 inpf.write(" Beta_n_Factor: %4.6e\n" % (beta_n_factor))
84
85def write_material_zone_part(inpf, zone_string, horizon, rho, K, G, Gc):
86 inpf.write(" Zone_%s:\n" % (zone_string))
87 inpf.write(" Type: PDState\n")
88 inpf.write(" Horizon: %4.6e\n" % (horizon))
89 inpf.write(" Density: %4.6e\n" % (rho))
90 inpf.write(" Compute_From_Classical: true\n")
91 inpf.write(" K: %4.6e\n" % (K))
92 inpf.write(" G: %4.6e\n" % (G))
93 inpf.write(" Gc: %4.6e\n" % (Gc))
94 inpf.write(" Influence_Function:\n")
95 inpf.write(" Type: 1\n")
96
97def copy_contact_zone(inpf, zone_numbers, zone_copy_from):
98 for s in zone_numbers:
99 inpf.write(" Zone_%d:\n" % (s))
100 inpf.write(" Copy_Contact_Data: " + print_int_list(zone_copy_from))
101
102def get_E(K, nu):
103 return 3. * K * (1. - 2. * nu)
104
105def get_G(E, nu):
106 return E / (2. * (1. + nu))
107
108
109def get_eff_k(k1, k2):
110 return 2. * k1 * k2 / (k1 + k2)
111
112
113def get_max(l):
114 l = np.array(l)
115 return l.max()
116
117def get_center(p1, p2):
118 return [0.5*(p1[i] + p2[i]) for i in range(3)]
119
120def rotate(p, theta, axis):
121
122 p_np = np.array(p)
123 axis_np = np.array(axis)
124
125 ct = np.cos(theta);
126 st = np.sin(theta);
127
128 # dot
129 p_dot_n = np.dot(p_np,axis_np)
130
131 # cross
132 n_cross_p = np.cross(axis_np, p_np)
133
134 return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p
135
136def get_ref_rect_points(center, L, W, add_center = False):
137
138 points = []
139 if add_center:
140 points.append(center)
141
142 points.append([center[0] - 0.5*L, center[1] - 0.5*W, center[2]])
143 points.append([center[0] + 0.5*L, center[1] - 0.5*W, center[2]])
144 points.append([center[0] + 0.5*L, center[1] + 0.5*W, center[2]])
145 points.append([center[0] - 0.5*L, center[1] + 0.5*W, center[2]])
146
147 return points
148
149def get_ref_triangle_points(center, radius, add_center = False):
150
151 # triangle
152 # 2
153 # +
154 #
155 #
156 # o
157 # 1
158 #
159 # +-------------------+
160 # 3 4
161
162 # center and radius
163 sim_Cx = center[0]
164 sim_Cy = center[1]
165 sim_Cz = center[2]
166 sim_radius = radius
167
168 cp = np.cos(np.pi/6.)
169 sp = np.sin(np.pi/6.)
170
171 points = []
172 if add_center:
173 points.append([sim_Cx, sim_Cy, sim_Cz])
174 points.append([sim_Cx, sim_Cy + sim_radius, sim_Cz])
175 points.append([sim_Cx - sim_radius*cp, sim_Cy - sim_radius*sp, sim_Cz])
176 points.append([sim_Cx + sim_radius*cp, sim_Cy - sim_radius*sp, sim_Cz])
177
178 return points
179
180def get_ref_hex_points(center, radius, add_center = False):
181
182 # drum2d
183 #
184 # v3 v2
185 # + +
186 #
187 #
188 # + o +
189 # v4 x v1
190 #
191 # + +
192 # v5 v6
193 #
194 # Axis is a vector from x to v1
195 #
196
197 axis = [1., 0., 0.]
198
199 # center and radius
200 sim_Cx = center[0]
201 sim_Cy = center[1]
202 sim_Cz = center[2]
203 sim_radius = radius
204
205 rotate_axis = [0., 0., 1.]
206
207 points = []
208 if add_center:
209 points.append(center)
210
211 for i in range(6):
212 xi = rotate(axis, i*np.pi/3., rotate_axis)
213 points.append([center[i] + radius * xi[i] for i in range(3)])
214
215 return points
216
217def get_ref_drum_points(center, radius, width, add_center = False):
218
219 # drum2d
220 #
221 # v3 v2
222 # + +
223 #
224 #
225 # + o +
226 # v4 x v1
227 #
228 # + +
229 # v5 v6
230 #
231 # Axis is a vector from x to v1
232 #
233
234 axis = [1., 0., 0.]
235
236 # center and radius
237 sim_Cx = center[0]
238 sim_Cy = center[1]
239 sim_Cz = center[2]
240 sim_radius = radius
241
242 rotate_axis = [0., 0., 1.]
243
244 x1 = rotate(axis, np.pi/3., rotate_axis)
245 x2 = rotate(axis, -np.pi/3., rotate_axis)
246
247 points = []
248 if add_center:
249 points.append(center)
250
251 # v1
252 points.append([center[i] + width*0.5*axis[i] for i in range(3)])
253 # v2
254 points.append([center[i] + radius*x1[i] for i in range(3)])
255 # v3
256 points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)])
257 # v4
258 points.append([center[i] - width*0.5*axis[i] for i in range(3)])
259 # v5
260 v6 = [center[i] + radius*x2[i] for i in range(3)]
261 points.append([v6[i] - radius*axis[i] for i in range(3)])
262 # v6
263 points.append(v6)
264
265 return points
266
267def does_rect_intersect_rect(r1, r2, padding):
268
269 # enlarge rectangle by adding padding
270 r1_padding = [r1[0] - padding, r1[1] - padding, r1[2], r1[3] + padding, r1[4] + padding, r1[5]]
271
272 return r1_padding[0] < r2[3] and r1_padding[3] > r2[0] and r1_padding[1] < r2[4] and r1_padding[4] > r2[1]
273
274def does_rect_intersect_rect_use_pair_coord(r1, r2, padding):
275
276 # enlarge rectangle by adding padding
277 p1_r1_padding = [r1[0][0] - padding, r1[0][1] - padding, r1[0][2]]
278 p2_r1_padding = [r1[1][0] + padding, r1[1][1] + padding, r1[0][2]]
279
280 p1_r2 = r2[0]
281 p2_r2 = r2[1]
282
283 return p1_r1_padding[0] < p2_r2[0] and p2_r1_padding[0] > p1_r2[0] and p1_r1_padding[1] < p2_r2[1] and p2_r1_padding[1] > p1_r2[1]
284
285
286def does_particle_intersect_rect(p, r2, padding):
287
288 pr = p[4]
289 pc = [p[1], p[2], p[3]]
290 p_rect = [pc[0] - pr, pc[1] - pr, pc[2], pc[1] + pr, pc[2] + pr, pc[2]]
291
292 return does_rect_intersect_rect(p_rect, r2, padding)
293
294def does_particle_intersect(p, particles, rect, padding):
295
296 # p -> [id, x, y, z, r]
297
298 p_center = [p[i+1] for i in range(3)]
299 p_r = p[4]
300 p_rect = [p_center[0] - p_r, p_center[1] - p_r, p_center[2], p_center[0] + p_r, p_center[1] + p_r, p_center[2]]
301
302 if p_rect[0] < rect[0] + padding or p_rect[1] < rect[1] + padding or p_rect[3] > rect[3] - padding or p_rect[4] > rect[4] - padding:
303 return True
304
305 for q in particles:
306 pq = np.array([p[i+1] - q[i+1] for i in range(3)])
307 if np.linalg.norm(pq) <= p[-1] + q[-1] + padding:
308 return True
309
310 return False
311
312
313def particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R2, id_choices1, id_choices2, N1, N2, rect, z_coord, add_orient = True):
314
315 np.random.seed(30)
316 sim_inp_dir = str(inp_dir)
317
318 """Generate particle location data"""
319
320 particles = []
321 points = []
322
323 method_to_use = 1
324
325 rect_L = rect[3] - rect[0]
326 rect_W = rect[4] - rect[1]
327
328 if method_to_use == 0:
329 pcount = 0
330 count = 0
331 select_count = 0
332 while pcount < N2 and count < 100*N2:
333 if count%N2 == 0:
334 print('large particle iter = ', count)
335
336 # random radius and center location
337 r = R2 + np.random.uniform(-0.1 * R2, 0.1 * R2)
338 x = center[0] + np.random.uniform(-0.5*rect_L + R2 + padding, 0.5*rect_L - R2- padding)
339 y = np.random.uniform(-0.5*rect_W + R2+padding, max_y - R2-padding)
340 p_zone = id_choices2[select_count % len(id_choices2)]
341 p = [p_zone, x, y, center[2], r]
342
343 # check if it collides of existing particles
344 pintersect = does_particle_intersect(p, particles, rect_L, rect_W, center, padding)
345
346 if pintersect == False:
347 particles.append(p)
348 pcount += 1
349 select_count += 1
350
351 count +=1
352
353 pcount = 0
354 count = 0
355 select_count = 0
356 while pcount < N1 and count < 100*N1:
357 if count%N1 == 0:
358 print('small particle iter = ', count)
359
360 # random radius and center location
361 r = R1 + np.random.uniform(-0.1 * R1, 0.1 * R1)
362 x = center[0] + np.random.uniform(-0.5*rect_L + R1 + padding, 0.5*rect_L - R1- padding)
363 y = np.random.uniform(-0.5*rect_W + R1 + padding, max_y - R1 - padding)
364 # p_zone = np.random.choice(id_choices1, size=1)[0]
365 # if np.random.uniform(0., 1.) > 0.25:
366 # p_zone = np.random.choice(id_choices1, size=1)[0]
367 p_zone = id_choices1[select_count % len(id_choices1)]
368 p = [p_zone, x,y,center[2], r]
369
370 # check if it collides of existing particles
371 pintersect = does_particle_intersect(p, particles, rect_L, rect_W, center, padding)
372
373 if pintersect == False:
374 particles.append(p)
375 pcount += 1
376 select_count += 1
377
378 count +=1
379
380 elif method_to_use == 1:
381
382 # find how approximate number of rows and cols we can have
383 check_r = R1
384 if R1 > R2:
385 check_r = R2
386
387 rows = int((max_y - rect[1]) / (2. * check_r))
388 cols = int(rect_L / (2. * check_r))
389
390 print(rows, cols, N1, N2)
391
392 rads = [R1, R2]
393
394 counter = 0
395 counter1 = 0
396 counter2 = 0
397 x_old = rect[0]
398 x_old_right = rect[3]
399 y_old = rect[1]
400
401 cx = 0.
402 cy = 0.
403 cz = center[2]
404 r = 0.
405
406
407 cy_accptd = []
408 cy_accptd.append(y_old)
409
410 row_rads = []
411 row_rads.append(R1)
412 row_rads.append(R2)
413
414 for i in range(rows):
415 print('row = {}'.format(i), flush=True)
416
417 if i > 0:
418 y_old = get_max(cy_accptd) + get_max(row_rads)
419
420 #y_old = i*(2*get_max([R1, R2])) + rect[1]
421
422 # print(i, cy, y_old)
423
424 # reset
425 row_rads = []
426 row_rads.append(R1)
427 row_rads.append(R2)
428
429 if y_old + padding + get_max(row_rads) < max_y:
430 num_p_cols = 0
431 j = 0
432 while True:
433 if num_p_cols > cols - 1 or j > 100*(N1 + N2):
434 break
435
436 if counter1 >= N1 and counter2 >= N2:
437 break
438
439 if j == 0:
440 # for odd row, take x_old as first and for even take as last
441 x_old = rect[0]
442 x_old_right = rect[3]
443
444 # type of particle (type 1 or 2, e.g., small or large)
445 p_type = np.random.choice([0,1], size=1)[0]
446 if counter1 >= N1:
447 p_type = 1
448 if counter2 >= N2 :
449 p_type = 0
450
451 r0 = rads[p_type]
452
453 # zone this particle belongs to
454 p_zone = 0
455 if p_type == 0:
456 p_zone = np.random.choice(id_choices1, size=1)[0]
457 else:
458 p_zone = np.random.choice(id_choices2, size=1)[0]
459
460 # perturb radius
461 r = r0 + np.random.uniform(-0.1 * r0, 0.1 * r0)
462
463 # for even row, we start from right edge of rectange
464 if i%2 == 0:
465 # random horizontal/vertical perturbation
466 rph = np.random.uniform(-0.1 * r0, 0.05 * r0)
467 rpv = np.random.uniform(-0.05 * r0, 0.05 * r0)
468 cx0 = x_old_right - padding - r
469 cx = cx0 - rph
470 cy = y_old + padding + r + rpv
471 else:
472 # random horizontal/vertical perturbation
473 rph = np.random.uniform(-0.05 * r0, 0.1 * r0)
474 rpv = np.random.uniform(-0.05 * r0, 0.05 * r0)
475 cx0 = x_old + padding + r
476 cx = cx0 + rph
477 cy = y_old + padding + r + rpv
478
479 # particle
480 p = [p_zone, cx, cy, cz, r]
481
482 # check for intersection
483 # if i < 3:
484 # print(p, particles, rect_L, rect_W, center, padding)
485 pintersect = does_particle_intersect(p, particles, rect, padding)
486
487 if pintersect == False:
488 print(i, p, flush=True)
489 particles.append(p)
490 row_rads.append(p[4])
491 cy_accptd.append(cy)
492 # set x_old
493 if i % 2 == 0:
494 x_old_right = cx - p[4]
495 else:
496 x_old = cx + p[4]
497
498 if p_type == 0:
499 counter1 += 1
500 elif p_type == 1:
501 counter2 += 1
502
503 counter += 1
504
505 num_p_cols += 1
506
507 j += 1
508
509 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
510 if add_orient:
511 inpf.write("i, x, y, z, r, o\n")
512 counter = 0
513 for p in particles:
514 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)))
515 counter += 1
516 else:
517 inpf.write("i, x, y, z, r\n")
518 for p in particles:
519 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
520
521
522 inpf.close()
523
524 # to visualize in paraview
525 write_vtu = False
526 if write_vtu:
527 points = []
528 rads = []
529 zones = []
530 for p in particles:
531 points.append([p[1], p[2], p[3]])
532 rads.append(p[-1])
533 zones.append(int(p[0]))
534
535 points = np.array(points)
536 rads = np.array(rads)
537 zones = np.array(zones)
538 mesh = pv.PolyData(points)
539 mesh["radius"] = rads
540 mesh["zone"] = zones
541 pv.save_meshio(sim_inp_dir + 'particle_locations_' + str(pp_tag) + '.vtu', mesh)
542
543 print('number of particles created = {}'.format(len(particles)), flush=True)
544
545 return len(particles), particles
546
547def particle_locations_old(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R2, id_choices1, id_choices2, N1, N2, rect, z_coord, add_orient = True):
548
549 """Generate particle location data"""
550
551
552 def does_intersect_old(p, r, R, particles, padding):
553
554 for q in particles:
555
556 pq = np.array([p[i] - q[i] for i in range(3)])
557 if np.linalg.norm(pq) <= r + R + padding:
558 return True
559
560 return False
561
562
563 def does_intersect_rect_old(p, r, particles, padding, rect):
564
565 # check intersection with rectangle
566 pr = [p[0] - r, p[1] - r, p[2], p[0] + r, p[1] + r, p[2]]
567
568 if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding:
569
570 # print('circle (xc = {}, r = {:5.3e}) intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
571
572 return True
573
574 # loop over particles
575 # for pi in particles:
576 # dx = [p[0] - pi[1], p[1] - pi[2], p[2] - pi[3]]
577 # rR = r + pi[4] + padding
578 # if np.linalg.norm(dx) < rR:
579 # return True
580
581 # print('circle (xc = {}, r = {:5.3e}) does not intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
582 return False
583
584 #
585 #
586 # Given rectangle, we want to fill it with particles of R1 and R2 radius
587 # in alternate fashion with padding between each particle
588 #
589 #
590 particles = []
591 points = []
592
593 # find how approximate number of rows and cols we can have
594 check_r = R1
595 if R1 > R2:
596 check_r = R2
597
598 rows = int((rect[3] - rect[0])/ check_r)
599 cols = int((rect[4] - rect[1])/ check_r)
600
601 rads = [R1, R2]
602
603 counter = 0
604 x_old = rect[0]
605 x_old_right = rect[3]
606 y_old = rect[1]
607 pad = padding
608
609 cx = 0.
610 cy = 0.
611 cz = 0.
612 r = 0.
613
614 row_rads = []
615 row_rads.append(R1)
616 row_rads.append(R2)
617
618 for i in range(rows):
619
620 if i > 0:
621 y_old = cy + get_max(row_rads)
622
623 # reset
624 row_rads = []
625 row_rads.append(R1)
626 row_rads.append(R2)
627
628 num_p_cols = 0
629 j = 0
630 while num_p_cols < cols - 1 and j < 500:
631
632 if j == 0:
633 # for odd row, take x_old as first and for even take as last
634 x_old = rect[0]
635 x_old_right = rect[3]
636
637 # type of particle (type 1 or 2, e.g., small or large)
638 p_type = counter % 2
639
640 r0 = rads[p_type]
641
642 # zone this particle belongs to
643 p_zone = 0
644 if p_type == 0:
645 p_zone = np.random.choice(id_choices1, size=1)[0]
646 else:
647 p_zone = np.random.choice(id_choices2, size=1)[0]
648
649 # perturb radius
650 r = r0 + np.random.uniform(-0.1 * r0, 0.1 * r0)
651 # r = r0
652 row_rads.append(r)
653
654 # random horizontal perturbation by 10% of radius
655 rph = np.random.uniform(-0.1 * r0, 0.1 * r0)
656
657 cx0 = x_old + pad + r
658 cx = cx0 + rph
659 # for even row, we start from right edge of rectange
660 if i%2 == 0:
661 cx0 = x_old_right - pad - r
662 cx = cx0 - rph
663
664 cy = y_old + pad + r
665 cz = rect[2]
666
667 inters = does_intersect_rect_old([cx, cy, cz], r, particles, pad, rect)
668 inters_parts = does_intersect_old
669
670 if inters == False:
671 particles.append([float(p_zone), cx, cy, cz, r])
672
673 # set x_old
674 if i % 2 == 0:
675 x_old_right = cx - r
676 else:
677 x_old = cx + r
678
679 counter += 1
680
681 num_p_cols += 1
682
683 j += 1
684
685
686 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
687 if add_orient:
688 inpf.write("i, x, y, z, r, o\n")
689 counter = 0
690 for p in particles:
691 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)))
692 counter += 1
693 else:
694 inpf.write("i, x, y, z, r\n")
695 for p in particles:
696 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
697
698
699 inpf.close()
700
701 # to visualize in paraview
702 write_vtu = False
703 if write_vtu:
704 points = []
705 rads = []
706 zones = []
707 for p in particles:
708 points.append([p[1], p[2], p[3]])
709 rads.append(p[-1])
710 zones.append(int(p[0]))
711
712 points = np.array(points)
713 rads = np.array(rads)
714 zones = np.array(zones)
715 mesh = pv.PolyData(points)
716 mesh["radius"] = rads
717 mesh["zone"] = zones
718 pv.save_meshio(sim_inp_dir + 'particle_locations_' + str(pp_tag) + '.vtu', mesh)
719
720 print('number of particles created = {}'.format(len(particles)))
721
722 return len(particles), particles
723
724
725def generate_cir_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
726
727 sim_inp_dir = str(inp_dir)
728
729 # center and radius
730 sim_Cx = center[0]
731 sim_Cy = center[1]
732 sim_Cz = center[2]
733 sim_radius = radius
734
735 # mesh size
736 sim_h = mesh_size
737
738 #
739 # create .geo file for gmsh
740 #
741 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
742 geof.write("cl__1 = 1;\n")
743 geof.write("Mesh.MshFileVersion = 2.2;\n")
744
745 #
746 # points
747 #
748 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz, sim_h));
749 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx + sim_radius, sim_Cy, sim_Cz, sim_h))
750 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx - sim_radius, sim_Cy, sim_Cz, sim_h))
751 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy + sim_radius, sim_Cz, sim_h))
752 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy - sim_radius, sim_Cz, sim_h))
753
754 #
755 # circlular arc
756 #
757 geof.write("Circle(1) = {2, 1, 4};\n")
758 geof.write("Circle(2) = {4, 1, 3};\n")
759 geof.write("Circle(3) = {3, 1, 5};\n")
760 geof.write("Circle(4) = {5, 1, 2};\n")
761
762 #
763 # surfaces
764 #
765 geof.write("Line Loop(1) = {2, 3, 4, 1};\n")
766
767 #
768 # plane surface
769 #
770 geof.write("Plane Surface(1) = {1};\n")
771
772 #
773 # physical surface
774 #
775 # tag = '"' + "a" + '"'
776 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
777
778 # # add center point to plane surface
779 geof.write("Point{1} In Surface {1};")
780
781 # close file
782 geof.close()
783
784
785def generate_hex_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
786
787 sim_inp_dir = str(inp_dir)
788
789 points = get_ref_hex_points(center, radius, True)
790
791 #
792 # create .geo file for gmsh
793 #
794 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
795 geof.write("cl__1 = 1;\n")
796 geof.write("Mesh.MshFileVersion = 2.2;\n")
797
798 #
799 # points
800 #
801 for i in range(7):
802 p = points[i]
803 sts = "Point({}) = ".format(i+1)
804 sts += "{"
805 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
806 sts += "};\n"
807 geof.write(sts);
808
809 #
810 # circlular arc
811 #
812 geof.write("Line(1) = {2, 3};\n")
813 geof.write("Line(2) = {3, 4};\n")
814 geof.write("Line(3) = {4, 5};\n")
815 geof.write("Line(4) = {5, 6};\n")
816 geof.write("Line(5) = {6, 7};\n")
817 geof.write("Line(6) = {7, 2};\n")
818
819 #
820 # surfaces
821 #
822 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
823
824 #
825 # plane surface
826 #
827 geof.write("Plane Surface(1) = {1};\n")
828
829 #
830 # physical surface
831 #
832 # tag = '"' + "a" + '"'
833 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
834
835 # # add center point to plane surface
836 geof.write("Point{1} In Surface {1};")
837
838 # close file
839 geof.close()
840
841
842def generate_tri_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
843
844 sim_inp_dir = str(inp_dir)
845
846 # points
847 points = get_ref_triangle_points(center, radius, True)
848
849 #
850 # create .geo file for gmsh
851 #
852 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
853 geof.write("cl__1 = 1;\n")
854 geof.write("Mesh.MshFileVersion = 2.2;\n")
855
856 #
857 # points
858 #
859 for i in range(4):
860 p = points[i]
861 sts = "Point({}) = ".format(i+1)
862 sts += "{"
863 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
864 sts += "};\n"
865 geof.write(sts);
866
867 #
868 # circlular arc
869 #
870 geof.write("Line(1) = {2, 3};\n")
871 geof.write("Line(2) = {3, 4};\n")
872 geof.write("Line(3) = {4, 2};\n")
873
874 #
875 # surfaces
876 #
877 geof.write("Line Loop(1) = {1, 2, 3};\n")
878
879 #
880 # plane surface
881 #
882 geof.write("Plane Surface(1) = {1};\n")
883
884 #
885 # physical surface
886 #
887 # tag = '"' + "a" + '"'
888 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
889
890 # # add center point to plane surface
891 geof.write("Point{1} In Surface {1};")
892
893 # close file
894 geof.close()
895
896
897def generate_drum2d_particle_gmsh_input(inp_dir, filename, center, radius, width, mesh_size, pp_tag):
898
899 sim_inp_dir = str(inp_dir)
900
901 # points
902 points = get_ref_drum_points(center, radius, width, True)
903
904 #
905 # create .geo file for gmsh
906 #
907 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
908 geof.write("cl__1 = 1;\n")
909 geof.write("Mesh.MshFileVersion = 2.2;\n")
910
911 #
912 # points
913 #
914 for i in range(7):
915 p = points[i]
916 sts = "Point({}) = ".format(i+1)
917 sts += "{"
918 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
919 sts += "};\n"
920 geof.write(sts);
921
922 #
923 # circlular arc
924 #
925 geof.write("Line(1) = {2, 3};\n")
926 geof.write("Line(2) = {3, 4};\n")
927 geof.write("Line(3) = {4, 5};\n")
928 geof.write("Line(4) = {5, 6};\n")
929 geof.write("Line(5) = {6, 7};\n")
930 geof.write("Line(6) = {7, 2};\n")
931
932 #
933 # surfaces
934 #
935 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
936
937 #
938 # plane surface
939 #
940 geof.write("Plane Surface(1) = {1};\n")
941
942 #
943 # physical surface
944 #
945 # tag = '"' + "a" + '"'
946 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
947
948 # # add center point to plane surface
949 geof.write("Point{1} In Surface {1};")
950
951 # close file
952 geof.close()
953
954def generate_rect_container_gmsh_input(inp_dir, filename, pi1, pi2, dx, dy, mesh_size, pp_tag):
955
956 sim_inp_dir = str(inp_dir)
957
958 innr_pts = get_ref_rect_points(get_center(pi1, pi2), pi2[0] - pi1[0], pi2[1] - pi1[1])
959
960 outr_pts = get_ref_rect_points(get_center(pi1, pi2), pi2[0] - pi1[0] + 2*dx, pi2[1] - pi1[1] + 2*dy)
961
962 h = mesh_size
963 hout = h # 5*h
964
965 #
966 # create .geo file for gmsh
967 #
968 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
969 geof.write("cl__1 = 1;\n")
970 geof.write("Mesh.MshFileVersion = 2.2;\n")
971
972 # points
973 p_id = 1
974 for i in range(4):
975 p_id = write_point_geo(geof, p_id, outr_pts[i], h)
976
977 for i in range(4):
978 p_id = write_point_geo(geof, p_id, innr_pts[i], h)
979
980 # line
981 l_id = 1
982 l_id = write_line_geo(geof, l_id, 1, 2)
983 l_id = write_line_geo(geof, l_id, 2, 3)
984 l_id = write_line_geo(geof, l_id, 3, 4)
985 l_id = write_line_geo(geof, l_id, 4, 1)
986
987 l_id = write_line_geo(geof, l_id, 5, 6)
988 l_id = write_line_geo(geof, l_id, 6, 7)
989 l_id = write_line_geo(geof, l_id, 7, 8)
990 l_id = write_line_geo(geof, l_id, 8, 5)
991
992 # line loop to define surface
993 geof.write("Line Loop(9) = {1, 2, 3, 4};\n")
994 geof.write("Line Loop(10) = {5, 6, 7, 8};\n")
995
996 # define surface
997 geof.write("Plane Surface(11) = {9, 10};\n")
998
999 # close file
1000 geof.close()
1001
1002def generate_rigid_rect_container_moving_wall_setup_gmsh_input(inp_dir, filename, outer_rect, inner_rect, mesh_size, pp_tag):
1003
1004 sim_inp_dir = str(inp_dir)
1005
1006 # outer rectangle
1007 sim_Lx_out1 = outer_rect[0]
1008 sim_Ly_out1 = outer_rect[1]
1009 sim_Lz_out1 = outer_rect[2]
1010 sim_Lx_out2 = outer_rect[3]
1011 sim_Ly_out2 = outer_rect[4]
1012 sim_Lz_out2 = outer_rect[5]
1013
1014 # inner rectangle
1015 sim_Lx_in1 = inner_rect[0]
1016 sim_Ly_in1 = inner_rect[1]
1017 sim_Lz_in1 = inner_rect[2]
1018 sim_Lx_in2 = inner_rect[3]
1019 sim_Ly_in2 = inner_rect[4]
1020 sim_Lz_in2 = inner_rect[5]
1021
1022 # mesh size
1023 sim_h = mesh_size
1024
1025 #
1026 # create .geo file for gmsh
1027 #
1028 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
1029 geof.write("cl__1 = 1;\n")
1030 geof.write("Mesh.MshFileVersion = 2.2;\n")
1031
1032 #
1033 # L7 L3
1034 # 4 + --- + 8 7 + --- + 3
1035 # | | | |
1036 # | | | |
1037 # | | L6 L4 | |
1038 # L8 | | | |
1039 # | | | | L2
1040 # | | L5 | |
1041 # | + ------------------- + |
1042 # | 5 6 |
1043 # | |
1044 # 1 + ------------------------------- + 2
1045 # L1
1046 #
1047 #
1048 # for point 8 and 7 --> choose y coordinate from outer rectangle and x
1049 # coordinate from inner rectangle
1050
1051 #
1052 # points
1053 #
1054 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h))
1055 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
1056 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
1057 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
1058
1059 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_in1, sim_Lz_in1, sim_h))
1060 geof.write("Point(6) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_in1, sim_Lz_in1, sim_h))
1061 geof.write("Point(7) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_out2, sim_Lz_in1, sim_h))
1062 geof.write("Point(8) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_out2, sim_Lz_in1, sim_h))
1063
1064 #
1065 # lines
1066 #
1067 geof.write("Line(1) = {1, 2};\n")
1068 geof.write("Line(2) = {2, 3};\n")
1069 geof.write("Line(3) = {3, 7};\n")
1070 geof.write("Line(4) = {7, 6};\n")
1071 geof.write("Line(5) = {6, 5};\n")
1072 geof.write("Line(6) = {5, 8};\n")
1073 geof.write("Line(7) = {8, 4};\n")
1074 geof.write("Line(8) = {4, 1};\n")
1075
1076 #
1077 # surfaces
1078 #
1079 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8};\n")
1080
1081 #
1082 # plane surface
1083 #
1084 geof.write("Plane Surface(1) = {1};\n")
1085
1086 #
1087 # physical surface
1088 #
1089 tag = '"' + "a" + '"'
1090 geof.write("Physical Surface(%s) = {1};\n" % (tag))
1091
1092 # close file
1093 geof.close()
1094
1095
1096def generate_moving_rect_wall_gmsh_input(inp_dir, filename, rectangle, mesh_size, pp_tag):
1097
1098 sim_inp_dir = str(inp_dir)
1099
1100 # outer rectangle
1101 sim_Lx_out1 = rectangle[0]
1102 sim_Ly_out1 = rectangle[1]
1103 sim_Lz_out1 = rectangle[2]
1104 sim_Lx_out2 = rectangle[3]
1105 sim_Ly_out2 = rectangle[4]
1106 sim_Lz_out2 = rectangle[5]
1107
1108 # mesh size
1109 sim_h = mesh_size
1110
1111 #
1112 # create .geo file for gmsh
1113 #
1114 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
1115 geof.write("cl__1 = 1;\n")
1116 geof.write("Mesh.MshFileVersion = 2.2;\n")
1117
1118 #
1119 # points
1120 #
1121 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h));
1122 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
1123 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
1124 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
1125
1126 #
1127 # lines
1128 #
1129 geof.write("Line(1) = {1, 2};\n")
1130 geof.write("Line(2) = {2, 3};\n")
1131 geof.write("Line(3) = {3, 4};\n")
1132 geof.write("Line(4) = {4, 1};\n")
1133
1134 #
1135 # surfaces
1136 #
1137 geof.write("Line Loop(1) = {1, 2, 3, 4};\n")
1138
1139 #
1140 # plane surface
1141 #
1142 geof.write("Plane Surface(1) = {1};\n")
1143
1144 #
1145 # physical surface
1146 #
1147 tag = '"' + "a" + '"'
1148 geof.write("Physical Surface(%s) = {1};\n" % (tag))
1149
1150 # close file
1151 geof.close()
1152
1153def create_input_file(inp_dir, pp_tag):
1154 """Generates input file for two-particle test"""
1155
1156 sim_inp_dir = str(inp_dir)
1157
1158 # 1 - small particle circle
1159 # 2 - small particle triangle
1160 # 3 - small particle drum2d
1161 # 4 - small particle hex
1162 # 5 - large particle circle
1163 # 6 - large particle triangle
1164 # 7 - large particle drum2d
1165 # 8 - large particle hex
1166 # 9 - rectangle wall rigid
1167 # 10 - top of wall that is moving
1168
1169
1170 center = [0., 0., 0.]
1171 R_small = 0.001
1172 R_large = 0.001
1173
1174 mesh_size = R_small / 5.
1175 horizon = 2 * mesh_size
1176
1177
1178 Lin, Win = 0.05, 0.04
1179 L, W = Lin + 1.5*horizon, Win + 1.5*horizon
1180
1181 in_rect = [center[0] - 0.5*Lin, center[1] - 0.5*Win, center[2], center[0] + 0.5*Lin, center[1] + 0.5*Win, center[2]]
1182 out_rect = [center[0] - 0.5*L, center[1] - 0.5*W, center[2], center[0] + 0.5*L, center[1] + 0.5*W, center[2]]
1183
1184 # moving wall
1185 moving_wall_y = 0.5*Win - 1.5*horizon
1186 moving_rect = [center[0] - 0.5*Lin, center[1] + moving_wall_y, center[2], center[0] + 0.5*Lin, center[1] + moving_wall_y + 1.5*horizon, center[2]]
1187
1188 # fixed container is annulus rectangle with top being removed
1189 # so we can represent it as outer rectangle minus rectangle that includes inner part and the top part
1190 remove_rect_from_out_rect = [in_rect[0], in_rect[1], in_rect[2], in_rect[3], out_rect[4], in_rect[5]]
1191 if moving_rect[4] > out_rect[4]:
1192 remove_rect_from_out_rect[4] = out_rect[4]
1193
1194 fixed_container_params = []
1195 for a in out_rect:
1196 fixed_container_params.append(a)
1197 for a in remove_rect_from_out_rect:
1198 fixed_container_params.append(a)
1199
1200 # container rectangle
1201 contain_rect = out_rect
1202
1203 # small circle
1204 small_circle = [R_small, center[0], center[1], center[2]]
1205 # small triangle
1206 small_triangle = small_circle
1207 # small drum2d
1208 w_small_drum2d = R_small * 0.2
1209 small_drum2d = [R_small, w_small_drum2d, center[0], center[1], center[2]]
1210 # small hex
1211 small_hex = small_circle
1212
1213 # large circle
1214 large_circle = [R_large, center[0], center[1], center[2]]
1215 # large triangle
1216 large_triangle = large_circle
1217 # large drum2d
1218 w_large_drum2d = R_large* 0.2
1219 large_drum2d = [R_large, w_large_drum2d, center[0], center[1], center[2]]
1220 # large hex
1221 large_hex = large_circle
1222
1223
1224
1225 final_time = 0.01
1226 num_steps = 4000
1227 # final_time = 0.00002
1228 # num_steps = 2000
1229 num_outputs = 10
1230 dt_out_n = num_steps / num_outputs
1231 test_dt_out_n = dt_out_n / 10
1232 perform_out = True
1233
1234
1235 rho_wall = 600.
1236 poisson_wall = 0.25
1237 K_wall = 1.e+4
1238 E_wall = get_E(K_wall, poisson_wall)
1239 G_wall = get_G(E_wall, poisson_wall)
1240 Gc_wall = 100.
1241
1242 rho_small = 600.
1243 poisson_small = poisson_wall
1244 K_small = 5.e+3
1245 E_small = get_E(K_small, poisson_small)
1246 G_small = get_G(E_small, poisson_small)
1247 Gc_small = 100.
1248
1249 rho_large = rho_small
1250 poisson_large = poisson_small
1251 K_large = K_small
1252 E_large = E_small
1253 G_large = G_small
1254 Gc_large = Gc_small
1255
1256
1259 R_contact_factor = 0.95
1260
1261 # Kn_V_max = 7.385158e+05
1262 # Kn = np.power(Kn_V_max, 2)
1263 # compute from bulk modulus
1264
1265 # from bulk modulus
1266 Kn_small_small = 18. * get_eff_k(K_small, K_small) / (np.pi * np.power(horizon, 5))
1267 Kn_large_large = 18. * get_eff_k(K_large, K_large) / (np.pi * np.power(horizon, 5))
1268 Kn_wall_wall = 18. * get_eff_k(K_wall, K_wall) / (np.pi * np.power(horizon, 5))
1269 Kn_small_large = 18. * get_eff_k(K_small, K_large) / (np.pi * np.power(horizon, 5))
1270 Kn_small_wall = 18. * get_eff_k(K_small, K_wall) / (np.pi * np.power(horizon, 5))
1271 Kn_large_wall = 18. * get_eff_k(K_large, K_wall) / (np.pi * np.power(horizon, 5))
1272
1273 # we do not want walls to interact
1274 Kn_wall_wall = 0.
1275
1276 beta_n_eps = 0.95
1277 friction_coeff = 0.5
1278 damping_active = False
1279 damping_active_floor = True
1280 friction_active = False
1281 beta_n_factor = 100.
1282 Kn_factor = 1.
1283
1284
1285 gravity_active = True
1286 gravity = [0., -10., 0.]
1287
1288
1289 neigh_search_factor = 10.
1290 neigh_search_interval = 100
1291 neigh_search_criteria = "simple_all"
1292
1293
1294 moving_wall_vert_velocity = -0.06
1295
1296
1299
1300 # generate particle locations
1301 padding = 1.1 * R_contact_factor * mesh_size
1302 max_y = moving_wall_y - 3*mesh_size
1303 # number of particles of small and large sizes
1304 N1, N2 = 300, 200
1305 id_choices1 = [0, 1, 2, 3]
1306 id_choices2 = [4, 5, 6, 7]
1307 num_particles_zones_1_to_8, particles_zones_1_to_8 = particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R_small, R_large, id_choices1, id_choices2, N1, N2, in_rect, z_coord = 0., add_orient = True)
1308
1309 # generate particle .geo file (large)
1310 zones_mesh_fnames = ["mesh_cir_small", "mesh_tri_small", "mesh_drum2d_small", "mesh_hex_small", "mesh_cir_large", "mesh_tri_large", "mesh_drum2d_large", "mesh_hex_large", "mesh_fixed_container", "mesh_moving_container"]
1311
1312
1313 generate_cir_particle_gmsh_input(inp_dir, zones_mesh_fnames[0], center, R_small, mesh_size, pp_tag)
1314 generate_cir_particle_gmsh_input(inp_dir, zones_mesh_fnames[4], center, R_large, mesh_size, pp_tag)
1315
1316
1317 generate_tri_particle_gmsh_input(inp_dir, zones_mesh_fnames[1], center, R_small, mesh_size, pp_tag)
1318 generate_tri_particle_gmsh_input(inp_dir, zones_mesh_fnames[5], center, R_large, mesh_size, pp_tag)
1319
1320
1321 generate_drum2d_particle_gmsh_input(inp_dir, zones_mesh_fnames[2], center, R_small, 2.*w_small_drum2d, mesh_size, pp_tag)
1322 generate_drum2d_particle_gmsh_input(inp_dir, zones_mesh_fnames[6], center, R_large, 2.*w_large_drum2d, mesh_size, pp_tag)
1323
1324
1325 generate_hex_particle_gmsh_input(inp_dir, zones_mesh_fnames[3], center, R_small, mesh_size, pp_tag)
1326 generate_hex_particle_gmsh_input(inp_dir, zones_mesh_fnames[7], center, R_large, mesh_size, pp_tag)
1327
1328
1329 generate_rigid_rect_container_moving_wall_setup_gmsh_input(inp_dir, zones_mesh_fnames[8], out_rect, in_rect, mesh_size, pp_tag)
1330
1331 generate_moving_rect_wall_gmsh_input(inp_dir, zones_mesh_fnames[9], moving_rect, mesh_size, pp_tag)
1332
1333 os.system("mkdir -p ../out")
1334
1335 for s in zones_mesh_fnames:
1336 print('\n')
1337 print(s)
1338 print("gmsh {}_{}.geo -2 &> /dev/null".format(s, pp_tag))
1339 print('\n')
1340
1341 os.system("gmsh {}_{}.geo -2".format(s, pp_tag))
1342 # os.system("gmsh {}_{}.geo -2 &> /dev/null".format(s, pp_tag))
1343 # os.system("gmsh {}_{}.geo -2 -o {}_{}.vtk &> /dev/null".format(s, pp_tag, s, pp_tag))
1344
1345
1346
1350 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
1351 inpf.write("Model:\n")
1352 inpf.write(" Dimension: 2\n")
1353 inpf.write(" Discretization_Type:\n")
1354 inpf.write(" Spatial: finite_difference\n")
1355 inpf.write(" Time: central_difference\n")
1356 inpf.write(" Final_Time: %4.6e\n" % (final_time))
1357 inpf.write(" Time_Steps: %d\n" % (num_steps))
1358
1359 # container info
1360 inpf.write("Container:\n")
1361 inpf.write(" Geometry:\n")
1362 inpf.write(" Type: rectangle\n")
1363 inpf.write(" Parameters: " + print_dbl_list(contain_rect))
1364
1365 # zone info
1366 inpf.write("Zone:\n")
1367 inpf.write(" Zones: 10\n")
1368
1369 for i in range(10):
1370 inpf.write(" Zone_%d:\n" % (i+1))
1371 if i > 7:
1372 inpf.write(" Is_Wall: true\n")
1373 else:
1374 inpf.write(" Is_Wall: false\n")
1375
1376 # particle info
1377 inpf.write("Particle:\n")
1378 inpf.write(" Test_Name: multi_particle_compressive_test\n")
1379
1380 particle_data = [['circle', small_circle], ['triangle', small_triangle], ['drum2d', small_drum2d], ['hexagon', small_hex], ['circle', large_circle], ['triangle', large_triangle], ['drum2d', large_drum2d], ['hexagon', large_hex]]
1381 for i in range(len(particle_data)):
1382 inpf.write(" Zone_%d:\n" % (i+1))
1383 inpf.write(" Type: %s\n" % (particle_data[i][0]))
1384 inpf.write(" Parameters: " + print_dbl_list((particle_data[i][1])))
1385
1386
1387 inpf.write(" Zone_9:\n")
1388 inpf.write(" Is_Wall: true\n")
1389 inpf.write(" Type: rectangle_minus_rectangle\n")
1390 inpf.write(" Parameters: " + print_dbl_list(fixed_container_params))
1391 inpf.write(" All_Dofs_Constrained: true\n")
1392 inpf.write(" Create_Particle_Using_ParticleZone_GeomObject: true\n")
1393
1394
1395 inpf.write(" Zone_10:\n")
1396 inpf.write(" Is_Wall: true\n")
1397 inpf.write(" Type: rectangle\n")
1398 inpf.write(" Parameters: " + print_dbl_list(moving_rect))
1399 inpf.write(" All_Dofs_Constrained: true\n")
1400 inpf.write(" Create_Particle_Using_ParticleZone_GeomObject: true\n")
1401
1402 # particle generation
1403 inpf.write("Particle_Generation:\n")
1404 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
1405 inpf.write(" File_Data_Type: loc_rad_orient\n")
1406
1407 # Mesh info
1408 inpf.write("Mesh:\n")
1409
1410 for i in range(10):
1411 inpf.write(" Zone_%d:\n" % (i+1))
1412 inpf.write(" File: %s\n" % (zones_mesh_fnames[i] + "_" + str(pp_tag) + ".msh"))
1413
1414 # Contact info
1415 inpf.write("Contact:\n")
1416
1417
1418 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "11", Kn_small_small)
1419
1420
1421 copy_contact_zone(inpf, [12, 13, 14], [1, 1])
1422
1423
1424 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "15", Kn_small_large)
1425
1426
1427 copy_contact_zone(inpf, [16, 17, 18], [1, 5])
1428
1429
1430 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "19", Kn_small_wall)
1431
1432
1433 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "110", Kn_small_wall)
1434
1435
1436 copy_contact_zone(inpf, [22, 23, 24], [1, 1])
1437
1438
1439 copy_contact_zone(inpf, [25, 26, 27, 28], [1, 5])
1440
1441
1442 copy_contact_zone(inpf, [29], [1, 9])
1443
1444
1445 copy_contact_zone(inpf, [210], [1, 10])
1446
1447
1448 copy_contact_zone(inpf, [33, 34], [1, 1])
1449
1450
1451 copy_contact_zone(inpf, [35, 36, 37, 38], [1, 5])
1452
1453
1454 copy_contact_zone(inpf, [39], [1, 9])
1455
1456
1457 copy_contact_zone(inpf, [310], [1, 10])
1458
1459
1460 copy_contact_zone(inpf, [44], [1, 1])
1461
1462
1463 copy_contact_zone(inpf, [45, 46, 47, 48], [1, 5])
1464
1465
1466 copy_contact_zone(inpf, [49], [1, 9])
1467
1468
1469 copy_contact_zone(inpf, [410], [1, 10])
1470
1471
1472 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "55", Kn_large_large)
1473
1474
1475 copy_contact_zone(inpf, [56, 57, 58], [5, 5])
1476
1477
1478 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "59", Kn_large_wall)
1479
1480
1481 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "510", Kn_large_wall)
1482
1483
1484 copy_contact_zone(inpf, [66, 67, 68], [5, 5])
1485
1486
1487 copy_contact_zone(inpf, [69], [5, 9])
1488
1489
1490 copy_contact_zone(inpf, [610], [5, 10])
1491
1492
1493 copy_contact_zone(inpf, [77, 78], [5, 5])
1494
1495
1496 copy_contact_zone(inpf, [79], [5, 9])
1497
1498
1499 copy_contact_zone(inpf, [710], [5, 10])
1500
1501
1502 copy_contact_zone(inpf, [88], [5, 5])
1503
1504
1505 copy_contact_zone(inpf, [89], [5, 9])
1506
1507
1508 copy_contact_zone(inpf, [810], [5, 10])
1509
1510
1511 write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "99", Kn_wall_wall)
1512
1513
1514 copy_contact_zone(inpf, [910], [9, 9])
1515
1516
1517 copy_contact_zone(inpf, [1010], [9, 9])
1518
1519 # Neighbor info
1520 inpf.write("Neighbor:\n")
1521 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
1522 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
1523 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
1524
1525 # Material info
1526 inpf.write("Material:\n")
1527
1528
1529 write_material_zone_part(inpf, "1", horizon, rho_small, K_small, G_small, Gc_small)
1530
1531
1532 inpf.write(" Zone_2:\n")
1533 inpf.write(" Copy_Material_Data: 1\n")
1534
1535
1536 inpf.write(" Zone_3:\n")
1537 inpf.write(" Copy_Material_Data: 1\n")
1538
1539
1540 inpf.write(" Zone_4:\n")
1541 inpf.write(" Copy_Material_Data: 1\n")
1542
1543
1544 write_material_zone_part(inpf, "5", horizon, rho_large, K_large, G_large, Gc_large)
1545
1546
1547 inpf.write(" Zone_6:\n")
1548 inpf.write(" Copy_Material_Data: 5\n")
1549
1550
1551 inpf.write(" Zone_7:\n")
1552 inpf.write(" Copy_Material_Data: 5\n")
1553
1554
1555 inpf.write(" Zone_8:\n")
1556 inpf.write(" Copy_Material_Data: 5\n")
1557
1558
1559 write_material_zone_part(inpf, "9", horizon, rho_wall, K_wall, G_wall, Gc_wall)
1560
1561
1562 inpf.write(" Zone_10:\n")
1563 inpf.write(" Copy_Material_Data: 9\n")
1564
1565 # Force
1566 if gravity_active == True:
1567 inpf.write("Force_BC:\n")
1568 inpf.write(" Gravity: " + print_dbl_list(gravity))
1569
1570 # Displacement
1571 inpf.write("Displacement_BC:\n")
1572 inpf.write(" Sets: 2\n")
1573
1574 inpf.write(" Set_1:\n")
1575 # wall particle id will be num_particles_zones_1_to_8
1576 inpf.write(" Particle_List: [%d]\n" % (num_particles_zones_1_to_8))
1577 inpf.write(" Direction: [1,2]\n")
1578 inpf.write(" Zero_Displacement: true\n")
1579
1580 inpf.write(" Set_2:\n")
1581 # wall particle id will be num_particles_zones_1_to_8 + 1
1582 inpf.write(" Particle_List: [%d]\n" % (num_particles_zones_1_to_8 + 1))
1583 inpf.write(" Direction: [2]\n")
1584 inpf.write(" Time_Function:\n")
1585 inpf.write(" Type: linear\n")
1586 inpf.write(" Parameters:\n")
1587 inpf.write(" - %4.6e\n" % (moving_wall_vert_velocity))
1588 inpf.write(" Spatial_Function:\n")
1589 inpf.write(" Type: constant\n")
1590
1591
1592 #
1593 # Output info
1594 #
1595 inpf.write("Output:\n")
1596 inpf.write(" Path: ../out/\n")
1597 inpf.write(" Tags:\n")
1598 inpf.write(" - Displacement\n")
1599 inpf.write(" - Velocity\n")
1600 inpf.write(" - Force\n")
1601 inpf.write(" - Force_Density\n")
1602 inpf.write(" - Damage_Z\n")
1603 inpf.write(" - Damage\n")
1604 inpf.write(" - Nodal_Volume\n")
1605 inpf.write(" - Zone_ID\n")
1606 inpf.write(" - Particle_ID\n")
1607 inpf.write(" - Fixity\n")
1608 inpf.write(" - Force_Fixity\n")
1609 inpf.write(" - Contact_Nodes\n")
1610 inpf.write(" - No_Fail_Node\n")
1611 inpf.write(" - Boundary_Node_Flag\n")
1612 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
1613 inpf.write(" Compress_Type: zlib\n")
1614 inpf.write(" Perform_FE_Out: false\n")
1615 if perform_out:
1616 inpf.write(" Perform_Out: true\n")
1617 else:
1618 inpf.write(" Perform_Out: false\n")
1619 inpf.write(" Test_Output_Interval: %d\n" % (test_dt_out_n))
1620
1621 inpf.write(" Debug: 2\n")
1622 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
1623
1624 # close file
1625 inpf.close()
1626
1627
1628
1630inp_dir = './'
1631pp_tag = 0
1632if len(sys.argv) > 1:
1633 pp_tag = int(sys.argv[1])
1634
1635create_input_file(inp_dir, pp_tag)