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 p_zone = id_choices1[counter1 % len(id_choices1)]
458 else:
459 # p_zone = np.random.choice(id_choices2, size=1)[0]
460 p_zone = id_choices2[counter2 % len(id_choices2)]
461
462 # perturb radius
463 r = r0 + np.random.uniform(-0.1 * r0, 0.1 * r0)
464
465 # for even row, we start from right edge of rectange
466 if i%2 == 0:
467 # random horizontal/vertical perturbation
468 rph = np.random.uniform(-0.1 * r0, 0.05 * r0)
469 rpv = np.random.uniform(-0.05 * r0, 0.05 * r0)
470 cx0 = x_old_right - padding - r
471 cx = cx0 - rph
472 cy = y_old + padding + r + rpv
473 else:
474 # random horizontal/vertical perturbation
475 rph = np.random.uniform(-0.05 * r0, 0.1 * r0)
476 rpv = np.random.uniform(-0.05 * r0, 0.05 * r0)
477 cx0 = x_old + padding + r
478 cx = cx0 + rph
479 cy = y_old + padding + r + rpv
480
481 # particle
482 p = [p_zone, cx, cy, cz, r]
483
484 # check for intersection
485 # if i < 3:
486 # print(p, particles, rect_L, rect_W, center, padding)
487 pintersect = does_particle_intersect(p, particles, rect, padding)
488
489 if pintersect == False:
490 print(i, p, flush=True)
491 particles.append(p)
492 row_rads.append(p[4])
493 cy_accptd.append(cy)
494 # set x_old
495 if i % 2 == 0:
496 x_old_right = cx - p[4]
497 else:
498 x_old = cx + p[4]
499
500 if p_type == 0:
501 counter1 += 1
502 elif p_type == 1:
503 counter2 += 1
504
505 counter += 1
506
507 num_p_cols += 1
508
509 j += 1
510
511 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
512 if add_orient:
513 inpf.write("i, x, y, z, r, o\n")
514 counter = 0
515 for p in particles:
516 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)))
517 counter += 1
518 else:
519 inpf.write("i, x, y, z, r\n")
520 for p in particles:
521 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
522
523
524 inpf.close()
525
526 # to visualize in paraview
527 write_vtu = False
528 if write_vtu:
529 points = []
530 rads = []
531 zones = []
532 for p in particles:
533 points.append([p[1], p[2], p[3]])
534 rads.append(p[-1])
535 zones.append(int(p[0]))
536
537 points = np.array(points)
538 rads = np.array(rads)
539 zones = np.array(zones)
540 mesh = pv.PolyData(points)
541 mesh["radius"] = rads
542 mesh["zone"] = zones
543 pv.save_meshio(sim_inp_dir + 'particle_locations_' + str(pp_tag) + '.vtu', mesh)
544
545 print('number of particles created = {}'.format(len(particles)), flush=True)
546
547 return len(particles), particles
548
549def 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):
550
551 """Generate particle location data"""
552
553
554 def does_intersect_old(p, r, R, particles, padding):
555
556 for q in particles:
557
558 pq = np.array([p[i] - q[i] for i in range(3)])
559 if np.linalg.norm(pq) <= r + R + padding:
560 return True
561
562 return False
563
564
565 def does_intersect_rect_old(p, r, particles, padding, rect):
566
567 # check intersection with rectangle
568 pr = [p[0] - r, p[1] - r, p[2], p[0] + r, p[1] + r, p[2]]
569
570 if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding:
571
572 # print('circle (xc = {}, r = {:5.3e}) intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
573
574 return True
575
576 # loop over particles
577 # for pi in particles:
578 # dx = [p[0] - pi[1], p[1] - pi[2], p[2] - pi[3]]
579 # rR = r + pi[4] + padding
580 # if np.linalg.norm(dx) < rR:
581 # return True
582
583 # print('circle (xc = {}, r = {:5.3e}) does not intersects rect = {} with pad = {:5.3e}'.format(p, r, rect, padding))
584 return False
585
586 #
587 #
588 # Given rectangle, we want to fill it with particles of R1 and R2 radius
589 # in alternate fashion with padding between each particle
590 #
591 #
592 particles = []
593 points = []
594
595 # find how approximate number of rows and cols we can have
596 check_r = R1
597 if R1 > R2:
598 check_r = R2
599
600 rows = int((rect[3] - rect[0])/ check_r)
601 cols = int((rect[4] - rect[1])/ check_r)
602
603 rads = [R1, R2]
604
605 counter = 0
606 x_old = rect[0]
607 x_old_right = rect[3]
608 y_old = rect[1]
609 pad = padding
610
611 cx = 0.
612 cy = 0.
613 cz = 0.
614 r = 0.
615
616 row_rads = []
617 row_rads.append(R1)
618 row_rads.append(R2)
619
620 for i in range(rows):
621
622 if i > 0:
623 y_old = cy + get_max(row_rads)
624
625 # reset
626 row_rads = []
627 row_rads.append(R1)
628 row_rads.append(R2)
629
630 num_p_cols = 0
631 j = 0
632 while num_p_cols < cols - 1 and j < 500:
633
634 if j == 0:
635 # for odd row, take x_old as first and for even take as last
636 x_old = rect[0]
637 x_old_right = rect[3]
638
639 # type of particle (type 1 or 2, e.g., small or large)
640 p_type = counter % 2
641
642 r0 = rads[p_type]
643
644 # zone this particle belongs to
645 p_zone = 0
646 if p_type == 0:
647 p_zone = np.random.choice(id_choices1, size=1)[0]
648 else:
649 p_zone = np.random.choice(id_choices2, size=1)[0]
650
651 # perturb radius
652 r = r0 + np.random.uniform(-0.1 * r0, 0.1 * r0)
653 # r = r0
654 row_rads.append(r)
655
656 # random horizontal perturbation by 10% of radius
657 rph = np.random.uniform(-0.1 * r0, 0.1 * r0)
658
659 cx0 = x_old + pad + r
660 cx = cx0 + rph
661 # for even row, we start from right edge of rectange
662 if i%2 == 0:
663 cx0 = x_old_right - pad - r
664 cx = cx0 - rph
665
666 cy = y_old + pad + r
667 cz = rect[2]
668
669 inters = does_intersect_rect_old([cx, cy, cz], r, particles, pad, rect)
670 inters_parts = does_intersect_old
671
672 if inters == False:
673 particles.append([float(p_zone), cx, cy, cz, r])
674
675 # set x_old
676 if i % 2 == 0:
677 x_old_right = cx - r
678 else:
679 x_old = cx + r
680
681 counter += 1
682
683 num_p_cols += 1
684
685 j += 1
686
687
688 inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
689 if add_orient:
690 inpf.write("i, x, y, z, r, o\n")
691 counter = 0
692 for p in particles:
693 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)))
694 counter += 1
695 else:
696 inpf.write("i, x, y, z, r\n")
697 for p in particles:
698 inpf.write("%d, %Lf, %Lf, %Lf, %Lf\n" % (int(p[0]), p[1], p[2], p[3], p[4]))
699
700
701 inpf.close()
702
703 # to visualize in paraview
704 write_vtu = False
705 if write_vtu:
706 points = []
707 rads = []
708 zones = []
709 for p in particles:
710 points.append([p[1], p[2], p[3]])
711 rads.append(p[-1])
712 zones.append(int(p[0]))
713
714 points = np.array(points)
715 rads = np.array(rads)
716 zones = np.array(zones)
717 mesh = pv.PolyData(points)
718 mesh["radius"] = rads
719 mesh["zone"] = zones
720 pv.save_meshio(sim_inp_dir + 'particle_locations_' + str(pp_tag) + '.vtu', mesh)
721
722 print('number of particles created = {}'.format(len(particles)))
723
724 return len(particles), particles
725
726
727def generate_cir_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
728
729 sim_inp_dir = str(inp_dir)
730
731 # center and radius
732 sim_Cx = center[0]
733 sim_Cy = center[1]
734 sim_Cz = center[2]
735 sim_radius = radius
736
737 # mesh size
738 sim_h = mesh_size
739
740 #
741 # create .geo file for gmsh
742 #
743 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
744 geof.write("cl__1 = 1;\n")
745 geof.write("Mesh.MshFileVersion = 2.2;\n")
746
747 #
748 # points
749 #
750 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy, sim_Cz, sim_h));
751 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx + sim_radius, sim_Cy, sim_Cz, sim_h))
752 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx - sim_radius, sim_Cy, sim_Cz, sim_h))
753 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy + sim_radius, sim_Cz, sim_h))
754 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Cx, sim_Cy - sim_radius, sim_Cz, sim_h))
755
756 #
757 # circlular arc
758 #
759 geof.write("Circle(1) = {2, 1, 4};\n")
760 geof.write("Circle(2) = {4, 1, 3};\n")
761 geof.write("Circle(3) = {3, 1, 5};\n")
762 geof.write("Circle(4) = {5, 1, 2};\n")
763
764 #
765 # surfaces
766 #
767 geof.write("Line Loop(1) = {2, 3, 4, 1};\n")
768
769 #
770 # plane surface
771 #
772 geof.write("Plane Surface(1) = {1};\n")
773
774 #
775 # physical surface
776 #
777 # tag = '"' + "a" + '"'
778 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
779
780 # # add center point to plane surface
781 geof.write("Point{1} In Surface {1};")
782
783 # close file
784 geof.close()
785
786
787def generate_hex_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
788
789 sim_inp_dir = str(inp_dir)
790
791 points = get_ref_hex_points(center, radius, True)
792
793 #
794 # create .geo file for gmsh
795 #
796 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
797 geof.write("cl__1 = 1;\n")
798 geof.write("Mesh.MshFileVersion = 2.2;\n")
799
800 #
801 # points
802 #
803 for i in range(7):
804 p = points[i]
805 sts = "Point({}) = ".format(i+1)
806 sts += "{"
807 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
808 sts += "};\n"
809 geof.write(sts);
810
811 #
812 # circlular arc
813 #
814 geof.write("Line(1) = {2, 3};\n")
815 geof.write("Line(2) = {3, 4};\n")
816 geof.write("Line(3) = {4, 5};\n")
817 geof.write("Line(4) = {5, 6};\n")
818 geof.write("Line(5) = {6, 7};\n")
819 geof.write("Line(6) = {7, 2};\n")
820
821 #
822 # surfaces
823 #
824 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
825
826 #
827 # plane surface
828 #
829 geof.write("Plane Surface(1) = {1};\n")
830
831 #
832 # physical surface
833 #
834 # tag = '"' + "a" + '"'
835 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
836
837 # # add center point to plane surface
838 geof.write("Point{1} In Surface {1};")
839
840 # close file
841 geof.close()
842
843
844def generate_tri_particle_gmsh_input(inp_dir, filename, center, radius, mesh_size, pp_tag):
845
846 sim_inp_dir = str(inp_dir)
847
848 # points
849 points = get_ref_triangle_points(center, radius, True)
850
851 #
852 # create .geo file for gmsh
853 #
854 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
855 geof.write("cl__1 = 1;\n")
856 geof.write("Mesh.MshFileVersion = 2.2;\n")
857
858 #
859 # points
860 #
861 for i in range(4):
862 p = points[i]
863 sts = "Point({}) = ".format(i+1)
864 sts += "{"
865 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
866 sts += "};\n"
867 geof.write(sts);
868
869 #
870 # circlular arc
871 #
872 geof.write("Line(1) = {2, 3};\n")
873 geof.write("Line(2) = {3, 4};\n")
874 geof.write("Line(3) = {4, 2};\n")
875
876 #
877 # surfaces
878 #
879 geof.write("Line Loop(1) = {1, 2, 3};\n")
880
881 #
882 # plane surface
883 #
884 geof.write("Plane Surface(1) = {1};\n")
885
886 #
887 # physical surface
888 #
889 # tag = '"' + "a" + '"'
890 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
891
892 # # add center point to plane surface
893 geof.write("Point{1} In Surface {1};")
894
895 # close file
896 geof.close()
897
898
899def generate_drum2d_particle_gmsh_input(inp_dir, filename, center, radius, width, mesh_size, pp_tag):
900
901 sim_inp_dir = str(inp_dir)
902
903 # points
904 points = get_ref_drum_points(center, radius, width, True)
905
906 #
907 # create .geo file for gmsh
908 #
909 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
910 geof.write("cl__1 = 1;\n")
911 geof.write("Mesh.MshFileVersion = 2.2;\n")
912
913 #
914 # points
915 #
916 for i in range(7):
917 p = points[i]
918 sts = "Point({}) = ".format(i+1)
919 sts += "{"
920 sts += "{}, {}, {}, {}".format(p[0], p[1], p[2], mesh_size)
921 sts += "};\n"
922 geof.write(sts);
923
924 #
925 # circlular arc
926 #
927 geof.write("Line(1) = {2, 3};\n")
928 geof.write("Line(2) = {3, 4};\n")
929 geof.write("Line(3) = {4, 5};\n")
930 geof.write("Line(4) = {5, 6};\n")
931 geof.write("Line(5) = {6, 7};\n")
932 geof.write("Line(6) = {7, 2};\n")
933
934 #
935 # surfaces
936 #
937 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6};\n")
938
939 #
940 # plane surface
941 #
942 geof.write("Plane Surface(1) = {1};\n")
943
944 #
945 # physical surface
946 #
947 # tag = '"' + "a" + '"'
948 # geof.write("Physical Surface(%s) = {1};\n" % (tag))
949
950 # # add center point to plane surface
951 geof.write("Point{1} In Surface {1};")
952
953 # close file
954 geof.close()
955
956def generate_rect_container_gmsh_input(inp_dir, filename, pi1, pi2, dx, dy, mesh_size, pp_tag):
957
958 sim_inp_dir = str(inp_dir)
959
960 innr_pts = get_ref_rect_points(get_center(pi1, pi2), pi2[0] - pi1[0], pi2[1] - pi1[1])
961
962 outr_pts = get_ref_rect_points(get_center(pi1, pi2), pi2[0] - pi1[0] + 2*dx, pi2[1] - pi1[1] + 2*dy)
963
964 h = mesh_size
965 hout = h # 5*h
966
967 #
968 # create .geo file for gmsh
969 #
970 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
971 geof.write("cl__1 = 1;\n")
972 geof.write("Mesh.MshFileVersion = 2.2;\n")
973
974 # points
975 p_id = 1
976 for i in range(4):
977 p_id = write_point_geo(geof, p_id, outr_pts[i], h)
978
979 for i in range(4):
980 p_id = write_point_geo(geof, p_id, innr_pts[i], h)
981
982 # line
983 l_id = 1
984 l_id = write_line_geo(geof, l_id, 1, 2)
985 l_id = write_line_geo(geof, l_id, 2, 3)
986 l_id = write_line_geo(geof, l_id, 3, 4)
987 l_id = write_line_geo(geof, l_id, 4, 1)
988
989 l_id = write_line_geo(geof, l_id, 5, 6)
990 l_id = write_line_geo(geof, l_id, 6, 7)
991 l_id = write_line_geo(geof, l_id, 7, 8)
992 l_id = write_line_geo(geof, l_id, 8, 5)
993
994 # line loop to define surface
995 geof.write("Line Loop(9) = {1, 2, 3, 4};\n")
996 geof.write("Line Loop(10) = {5, 6, 7, 8};\n")
997
998 # define surface
999 geof.write("Plane Surface(11) = {9, 10};\n")
1000
1001 # close file
1002 geof.close()
1003
1004def generate_rigid_rect_container_moving_wall_setup_gmsh_input(inp_dir, filename, outer_rect, inner_rect, mesh_size, pp_tag):
1005
1006 sim_inp_dir = str(inp_dir)
1007
1008 # outer rectangle
1009 sim_Lx_out1 = outer_rect[0]
1010 sim_Ly_out1 = outer_rect[1]
1011 sim_Lz_out1 = outer_rect[2]
1012 sim_Lx_out2 = outer_rect[3]
1013 sim_Ly_out2 = outer_rect[4]
1014 sim_Lz_out2 = outer_rect[5]
1015
1016 # inner rectangle
1017 sim_Lx_in1 = inner_rect[0]
1018 sim_Ly_in1 = inner_rect[1]
1019 sim_Lz_in1 = inner_rect[2]
1020 sim_Lx_in2 = inner_rect[3]
1021 sim_Ly_in2 = inner_rect[4]
1022 sim_Lz_in2 = inner_rect[5]
1023
1024 # mesh size
1025 sim_h = mesh_size
1026
1027 #
1028 # create .geo file for gmsh
1029 #
1030 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
1031 geof.write("cl__1 = 1;\n")
1032 geof.write("Mesh.MshFileVersion = 2.2;\n")
1033
1034 #
1035 # L7 L3
1036 # 4 + --- + 8 7 + --- + 3
1037 # | | | |
1038 # | | | |
1039 # | | L6 L4 | |
1040 # L8 | | | |
1041 # | | | | L2
1042 # | | L5 | |
1043 # | + ------------------- + |
1044 # | 5 6 |
1045 # | |
1046 # 1 + ------------------------------- + 2
1047 # L1
1048 #
1049 #
1050 # for point 8 and 7 --> choose y coordinate from outer rectangle and x
1051 # coordinate from inner rectangle
1052
1053 #
1054 # points
1055 #
1056 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h))
1057 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
1058 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
1059 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
1060
1061 geof.write("Point(5) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_in1, sim_Lz_in1, sim_h))
1062 geof.write("Point(6) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_in1, sim_Lz_in1, sim_h))
1063 geof.write("Point(7) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in2, sim_Ly_out2, sim_Lz_in1, sim_h))
1064 geof.write("Point(8) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_in1, sim_Ly_out2, sim_Lz_in1, sim_h))
1065
1066 #
1067 # lines
1068 #
1069 geof.write("Line(1) = {1, 2};\n")
1070 geof.write("Line(2) = {2, 3};\n")
1071 geof.write("Line(3) = {3, 7};\n")
1072 geof.write("Line(4) = {7, 6};\n")
1073 geof.write("Line(5) = {6, 5};\n")
1074 geof.write("Line(6) = {5, 8};\n")
1075 geof.write("Line(7) = {8, 4};\n")
1076 geof.write("Line(8) = {4, 1};\n")
1077
1078 #
1079 # surfaces
1080 #
1081 geof.write("Line Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8};\n")
1082
1083 #
1084 # plane surface
1085 #
1086 geof.write("Plane Surface(1) = {1};\n")
1087
1088 #
1089 # physical surface
1090 #
1091 tag = '"' + "a" + '"'
1092 geof.write("Physical Surface(%s) = {1};\n" % (tag))
1093
1094 # close file
1095 geof.close()
1096
1097
1098def generate_moving_rect_wall_gmsh_input(inp_dir, filename, rectangle, mesh_size, pp_tag):
1099
1100 sim_inp_dir = str(inp_dir)
1101
1102 # outer rectangle
1103 sim_Lx_out1 = rectangle[0]
1104 sim_Ly_out1 = rectangle[1]
1105 sim_Lz_out1 = rectangle[2]
1106 sim_Lx_out2 = rectangle[3]
1107 sim_Ly_out2 = rectangle[4]
1108 sim_Lz_out2 = rectangle[5]
1109
1110 # mesh size
1111 sim_h = mesh_size
1112
1113 #
1114 # create .geo file for gmsh
1115 #
1116 geof = open(sim_inp_dir + filename + '_' + str(pp_tag) + '.geo','w')
1117 geof.write("cl__1 = 1;\n")
1118 geof.write("Mesh.MshFileVersion = 2.2;\n")
1119
1120 #
1121 # points
1122 #
1123 geof.write("Point(1) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out1, sim_Lz_out1, sim_h));
1124 geof.write("Point(2) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out1, sim_Lz_out1, sim_h))
1125 geof.write("Point(3) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out2, sim_Ly_out2, sim_Lz_out1, sim_h))
1126 geof.write("Point(4) = {%4.6e, %4.6e, %4.6e, %4.6e};\n" % (sim_Lx_out1, sim_Ly_out2, sim_Lz_out1, sim_h))
1127
1128 #
1129 # lines
1130 #
1131 geof.write("Line(1) = {1, 2};\n")
1132 geof.write("Line(2) = {2, 3};\n")
1133 geof.write("Line(3) = {3, 4};\n")
1134 geof.write("Line(4) = {4, 1};\n")
1135
1136 #
1137 # surfaces
1138 #
1139 geof.write("Line Loop(1) = {1, 2, 3, 4};\n")
1140
1141 #
1142 # plane surface
1143 #
1144 geof.write("Plane Surface(1) = {1};\n")
1145
1146 #
1147 # physical surface
1148 #
1149 tag = '"' + "a" + '"'
1150 geof.write("Physical Surface(%s) = {1};\n" % (tag))
1151
1152 # close file
1153 geof.close()
1154
1155def create_input_file(inp_dir, pp_tag):
1156 """Generates input file for two-particle test"""
1157
1158 sim_inp_dir = str(inp_dir)
1159
1160 # 1 - small particle circle
1161 # 2 - small particle triangle
1162 # 3 - small particle drum2d
1163 # 4 - small particle hex
1164 # 5 - large particle circle
1165 # 6 - large particle triangle
1166 # 7 - large particle drum2d
1167 # 8 - large particle hex
1168 # 9 - rectangle wall rigid
1169 # 10 - top of wall that is moving
1170
1171
1172 center = [0., 0., 0.]
1173 R_small = 0.001
1174 R_large = 0.001
1175
1176 mesh_size = R_small / 5.
1177 horizon = 2 * mesh_size
1178
1179
1180 Lin, Win = 0.012, 0.008
1181 L, W = Lin + 1.5*horizon, Win + 1.5*horizon
1182
1183 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]]
1184 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]]
1185
1186 # moving wall
1187 moving_wall_y = 0.5*Win - 1.5*horizon
1188 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]]
1189
1190 # fixed container is annulus rectangle with top being removed
1191 # so we can represent it as outer rectangle minus rectangle that includes inner part and the top part
1192 remove_rect_from_out_rect = [in_rect[0], in_rect[1], in_rect[2], in_rect[3], out_rect[4], in_rect[5]]
1193 if moving_rect[4] > out_rect[4]:
1194 remove_rect_from_out_rect[4] = out_rect[4]
1195
1196 fixed_container_params = []
1197 for a in out_rect:
1198 fixed_container_params.append(a)
1199 for a in remove_rect_from_out_rect:
1200 fixed_container_params.append(a)
1201
1202 # container rectangle
1203 contain_rect = out_rect
1204
1205 # small circle
1206 small_circle = [R_small, center[0], center[1], center[2]]
1207 # small triangle
1208 small_triangle = small_circle
1209 # small drum2d
1210 w_small_drum2d = R_small * 0.2
1211 small_drum2d = [R_small, w_small_drum2d, center[0], center[1], center[2]]
1212 # small hex
1213 small_hex = small_circle
1214
1215 # large circle
1216 large_circle = [R_large, center[0], center[1], center[2]]
1217 # large triangle
1218 large_triangle = large_circle
1219 # large drum2d
1220 w_large_drum2d = R_large* 0.2
1221 large_drum2d = [R_large, w_large_drum2d, center[0], center[1], center[2]]
1222 # large hex
1223 large_hex = large_circle
1224
1225
1226
1227 final_time = 0.05
1228 num_steps = 20000
1229 # final_time = 0.00002
1230 # num_steps = 2000
1231 num_outputs = 10
1232 dt_out_n = num_steps / num_outputs
1233 test_dt_out_n = dt_out_n / 10
1234 perform_out = True
1235
1236
1237 rho_wall = 600.
1238 poisson_wall = 0.25
1239 K_wall = 1.e+4
1240 E_wall = get_E(K_wall, poisson_wall)
1241 G_wall = get_G(E_wall, poisson_wall)
1242 Gc_wall = 100.
1243
1244 rho_small = 600.
1245 poisson_small = poisson_wall
1246 K_small = 5.e+3
1247 E_small = get_E(K_small, poisson_small)
1248 G_small = get_G(E_small, poisson_small)
1249 Gc_small = 100.
1250
1251 rho_large = rho_small
1252 poisson_large = poisson_small
1253 K_large = K_small
1254 E_large = E_small
1255 G_large = G_small
1256 Gc_large = Gc_small
1257
1258
1261 R_contact_factor = 0.95
1262
1263 # Kn_V_max = 7.385158e+05
1264 # Kn = np.power(Kn_V_max, 2)
1265 # compute from bulk modulus
1266
1267 # from bulk modulus
1268 Kn_small_small = 18. * get_eff_k(K_small, K_small) / (np.pi * np.power(horizon, 5))
1269 Kn_large_large = 18. * get_eff_k(K_large, K_large) / (np.pi * np.power(horizon, 5))
1270 Kn_wall_wall = 18. * get_eff_k(K_wall, K_wall) / (np.pi * np.power(horizon, 5))
1271 Kn_small_large = 18. * get_eff_k(K_small, K_large) / (np.pi * np.power(horizon, 5))
1272 Kn_small_wall = 18. * get_eff_k(K_small, K_wall) / (np.pi * np.power(horizon, 5))
1273 Kn_large_wall = 18. * get_eff_k(K_large, K_wall) / (np.pi * np.power(horizon, 5))
1274
1275 # we do not want walls to interact
1276 Kn_wall_wall = 0.
1277
1278 beta_n_eps = 0.95
1279 friction_coeff = 0.5
1280 damping_active = False
1281 damping_active_floor = True
1282 friction_active = False
1283 beta_n_factor = 100.
1284 Kn_factor = 1.
1285
1286
1287 gravity_active = True
1288 gravity = [0., -10., 0.]
1289
1290
1291 neigh_search_factor = 10.
1292 neigh_search_interval = 100
1293 neigh_search_criteria = "simple_all"
1294
1295
1296 moving_wall_vert_velocity = -0.06
1297
1298
1301
1302 # generate particle locations
1303 padding = 1.1 * R_contact_factor * mesh_size
1304 max_y = moving_wall_y - 3*mesh_size
1305 # number of particles of small and large sizes
1306 N1, N2 = 10, 8
1307 id_choices1 = [0, 1, 2, 3]
1308 id_choices2 = [4, 5, 6, 7]
1309 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)
1310
1311 # generate particle .geo file (large)
1312 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"]
1313
1314
1315 generate_cir_particle_gmsh_input(inp_dir, zones_mesh_fnames[0], center, R_small, mesh_size, pp_tag)
1316 generate_cir_particle_gmsh_input(inp_dir, zones_mesh_fnames[4], center, R_large, mesh_size, pp_tag)
1317
1318
1319 generate_tri_particle_gmsh_input(inp_dir, zones_mesh_fnames[1], center, R_small, mesh_size, pp_tag)
1320 generate_tri_particle_gmsh_input(inp_dir, zones_mesh_fnames[5], center, R_large, mesh_size, pp_tag)
1321
1322
1323 generate_drum2d_particle_gmsh_input(inp_dir, zones_mesh_fnames[2], center, R_small, 2.*w_small_drum2d, mesh_size, pp_tag)
1324 generate_drum2d_particle_gmsh_input(inp_dir, zones_mesh_fnames[6], center, R_large, 2.*w_large_drum2d, mesh_size, pp_tag)
1325
1326
1327 generate_hex_particle_gmsh_input(inp_dir, zones_mesh_fnames[3], center, R_small, mesh_size, pp_tag)
1328 generate_hex_particle_gmsh_input(inp_dir, zones_mesh_fnames[7], center, R_large, mesh_size, pp_tag)
1329
1330
1331 generate_rigid_rect_container_moving_wall_setup_gmsh_input(inp_dir, zones_mesh_fnames[8], out_rect, in_rect, mesh_size, pp_tag)
1332
1333 generate_moving_rect_wall_gmsh_input(inp_dir, zones_mesh_fnames[9], moving_rect, mesh_size, pp_tag)
1334
1335 os.system("mkdir -p ../out")
1336
1337 for s in zones_mesh_fnames:
1338 print('\n')
1339 print(s)
1340 print("gmsh {}_{}.geo -2 &> /dev/null".format(s, pp_tag))
1341 print('\n')
1342
1343 os.system("gmsh {}_{}.geo -2".format(s, pp_tag))
1344 # os.system("gmsh {}_{}.geo -2 &> /dev/null".format(s, pp_tag))
1345 # os.system("gmsh {}_{}.geo -2 -o {}_{}.vtk &> /dev/null".format(s, pp_tag, s, pp_tag))
1346
1347
1348
1352 inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
1353 inpf.write("Model:\n")
1354 inpf.write(" Dimension: 2\n")
1355 inpf.write(" Discretization_Type:\n")
1356 inpf.write(" Spatial: finite_difference\n")
1357 inpf.write(" Time: central_difference\n")
1358 inpf.write(" Final_Time: %4.6e\n" % (final_time))
1359 inpf.write(" Time_Steps: %d\n" % (num_steps))
1360
1361 # container info
1362 inpf.write("Container:\n")
1363 inpf.write(" Geometry:\n")
1364 inpf.write(" Type: rectangle\n")
1365 inpf.write(" Parameters: " + print_dbl_list(contain_rect))
1366
1367 # zone info
1368 inpf.write("Zone:\n")
1369 inpf.write(" Zones: 10\n")
1370
1371 for i in range(10):
1372 inpf.write(" Zone_%d:\n" % (i+1))
1373 if i > 7:
1374 inpf.write(" Is_Wall: true\n")
1375 else:
1376 inpf.write(" Is_Wall: false\n")
1377
1378 # particle info
1379 inpf.write("Particle:\n")
1380 inpf.write(" Test_Name: multi_particle_compressive_test\n")
1381
1382 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]]
1383 for i in range(len(particle_data)):
1384 inpf.write(" Zone_%d:\n" % (i+1))
1385 inpf.write(" Type: %s\n" % (particle_data[i][0]))
1386 inpf.write(" Parameters: " + print_dbl_list((particle_data[i][1])))
1387
1388
1389 inpf.write(" Zone_9:\n")
1390 inpf.write(" Is_Wall: true\n")
1391 inpf.write(" Type: rectangle_minus_rectangle\n")
1392 inpf.write(" Parameters: " + print_dbl_list(fixed_container_params))
1393 inpf.write(" All_Dofs_Constrained: true\n")
1394 inpf.write(" Create_Particle_Using_ParticleZone_GeomObject: true\n")
1395
1396
1397 inpf.write(" Zone_10:\n")
1398 inpf.write(" Is_Wall: true\n")
1399 inpf.write(" Type: rectangle\n")
1400 inpf.write(" Parameters: " + print_dbl_list(moving_rect))
1401 inpf.write(" All_Dofs_Constrained: true\n")
1402 inpf.write(" Create_Particle_Using_ParticleZone_GeomObject: true\n")
1403
1404 # particle generation
1405 inpf.write("Particle_Generation:\n")
1406 inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
1407 inpf.write(" File_Data_Type: loc_rad_orient\n")
1408
1409 # Mesh info
1410 inpf.write("Mesh:\n")
1411
1412 for i in range(10):
1413 inpf.write(" Zone_%d:\n" % (i+1))
1414 inpf.write(" File: %s\n" % (zones_mesh_fnames[i] + "_" + str(pp_tag) + ".msh"))
1415
1416 # Contact info
1417 inpf.write("Contact:\n")
1418
1419
1420 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)
1421
1422
1423 copy_contact_zone(inpf, [12, 13, 14], [1, 1])
1424
1425
1426 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)
1427
1428
1429 copy_contact_zone(inpf, [16, 17, 18], [1, 5])
1430
1431
1432 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)
1433
1434
1435 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)
1436
1437
1438 copy_contact_zone(inpf, [22, 23, 24], [1, 1])
1439
1440
1441 copy_contact_zone(inpf, [25, 26, 27, 28], [1, 5])
1442
1443
1444 copy_contact_zone(inpf, [29], [1, 9])
1445
1446
1447 copy_contact_zone(inpf, [210], [1, 10])
1448
1449
1450 copy_contact_zone(inpf, [33, 34], [1, 1])
1451
1452
1453 copy_contact_zone(inpf, [35, 36, 37, 38], [1, 5])
1454
1455
1456 copy_contact_zone(inpf, [39], [1, 9])
1457
1458
1459 copy_contact_zone(inpf, [310], [1, 10])
1460
1461
1462 copy_contact_zone(inpf, [44], [1, 1])
1463
1464
1465 copy_contact_zone(inpf, [45, 46, 47, 48], [1, 5])
1466
1467
1468 copy_contact_zone(inpf, [49], [1, 9])
1469
1470
1471 copy_contact_zone(inpf, [410], [1, 10])
1472
1473
1474 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)
1475
1476
1477 copy_contact_zone(inpf, [56, 57, 58], [5, 5])
1478
1479
1480 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)
1481
1482
1483 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)
1484
1485
1486 copy_contact_zone(inpf, [66, 67, 68], [5, 5])
1487
1488
1489 copy_contact_zone(inpf, [69], [5, 9])
1490
1491
1492 copy_contact_zone(inpf, [610], [5, 10])
1493
1494
1495 copy_contact_zone(inpf, [77, 78], [5, 5])
1496
1497
1498 copy_contact_zone(inpf, [79], [5, 9])
1499
1500
1501 copy_contact_zone(inpf, [710], [5, 10])
1502
1503
1504 copy_contact_zone(inpf, [88], [5, 5])
1505
1506
1507 copy_contact_zone(inpf, [89], [5, 9])
1508
1509
1510 copy_contact_zone(inpf, [810], [5, 10])
1511
1512
1513 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)
1514
1515
1516 copy_contact_zone(inpf, [910], [9, 9])
1517
1518
1519 copy_contact_zone(inpf, [1010], [9, 9])
1520
1521 # Neighbor info
1522 inpf.write("Neighbor:\n")
1523 inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
1524 inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
1525 inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))
1526
1527 # Material info
1528 inpf.write("Material:\n")
1529
1530
1531 write_material_zone_part(inpf, "1", horizon, rho_small, K_small, G_small, Gc_small)
1532
1533
1534 inpf.write(" Zone_2:\n")
1535 inpf.write(" Copy_Material_Data: 1\n")
1536
1537
1538 inpf.write(" Zone_3:\n")
1539 inpf.write(" Copy_Material_Data: 1\n")
1540
1541
1542 inpf.write(" Zone_4:\n")
1543 inpf.write(" Copy_Material_Data: 1\n")
1544
1545
1546 write_material_zone_part(inpf, "5", horizon, rho_large, K_large, G_large, Gc_large)
1547
1548
1549 inpf.write(" Zone_6:\n")
1550 inpf.write(" Copy_Material_Data: 5\n")
1551
1552
1553 inpf.write(" Zone_7:\n")
1554 inpf.write(" Copy_Material_Data: 5\n")
1555
1556
1557 inpf.write(" Zone_8:\n")
1558 inpf.write(" Copy_Material_Data: 5\n")
1559
1560
1561 write_material_zone_part(inpf, "9", horizon, rho_wall, K_wall, G_wall, Gc_wall)
1562
1563
1564 inpf.write(" Zone_10:\n")
1565 inpf.write(" Copy_Material_Data: 9\n")
1566
1567 # Force
1568 if gravity_active == True:
1569 inpf.write("Force_BC:\n")
1570 inpf.write(" Gravity: " + print_dbl_list(gravity))
1571
1572 # Displacement
1573 inpf.write("Displacement_BC:\n")
1574 inpf.write(" Sets: 2\n")
1575
1576 inpf.write(" Set_1:\n")
1577 # wall particle id will be num_particles_zones_1_to_8
1578 inpf.write(" Particle_List: [%d]\n" % (num_particles_zones_1_to_8))
1579 inpf.write(" Direction: [1,2]\n")
1580 inpf.write(" Zero_Displacement: true\n")
1581
1582 inpf.write(" Set_2:\n")
1583 # wall particle id will be num_particles_zones_1_to_8 + 1
1584 inpf.write(" Particle_List: [%d]\n" % (num_particles_zones_1_to_8 + 1))
1585 inpf.write(" Direction: [2]\n")
1586 inpf.write(" Time_Function:\n")
1587 inpf.write(" Type: linear\n")
1588 inpf.write(" Parameters:\n")
1589 inpf.write(" - %4.6e\n" % (moving_wall_vert_velocity))
1590 inpf.write(" Spatial_Function:\n")
1591 inpf.write(" Type: constant\n")
1592
1593
1594 #
1595 # Output info
1596 #
1597 inpf.write("Output:\n")
1598 inpf.write(" Path: ../out/\n")
1599 inpf.write(" Tags:\n")
1600 inpf.write(" - Displacement\n")
1601 inpf.write(" - Velocity\n")
1602 inpf.write(" - Force\n")
1603 inpf.write(" - Force_Density\n")
1604 inpf.write(" - Damage_Z\n")
1605 inpf.write(" - Damage\n")
1606 inpf.write(" - Nodal_Volume\n")
1607 inpf.write(" - Zone_ID\n")
1608 inpf.write(" - Particle_ID\n")
1609 inpf.write(" - Fixity\n")
1610 inpf.write(" - Force_Fixity\n")
1611 inpf.write(" - Contact_Nodes\n")
1612 inpf.write(" - No_Fail_Node\n")
1613 inpf.write(" - Boundary_Node_Flag\n")
1614 inpf.write(" Output_Interval: %d\n" % (dt_out_n))
1615 inpf.write(" Compress_Type: zlib\n")
1616 inpf.write(" Perform_FE_Out: false\n")
1617 if perform_out:
1618 inpf.write(" Perform_Out: true\n")
1619 else:
1620 inpf.write(" Perform_Out: false\n")
1621 inpf.write(" Test_Output_Interval: %d\n" % (test_dt_out_n))
1622
1623 inpf.write(" Debug: 3\n")
1624 inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))
1625
1626 # close file
1627 inpf.close()
1628
1629
1630
1632inp_dir = './'
1633pp_tag = 0
1634if len(sys.argv) > 1:
1635 pp_tag = int(sys.argv[1])
1636
1637create_input_file(inp_dir, pp_tag)