tpacking.c - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tpacking.c (2824B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <err.h>
4 #include <math.h>
5 #include "util.h"
6 #include "granular.h"
7 #include "grain.h"
8 #include "arrays.h"
9
10 size_t
11 rectangular_packing(struct grain **grains,
12 size_t n[3],
13 double diameter_min, double diameter_max,
14 double (*gsd)(double min, double max),
15 double padding_factor,
16 double origo[3])
17 {
18 size_t i, j, k, l, ng, N = 0;
19 double dx_padding = diameter_max * padding_factor;
20 double dx = diameter_max + dx_padding;
21
22 if (diameter_max < diameter_min)
23 errx(1, "%s: diameter_max (%g) is smaller than diameter_min (%g)",
24 __func__, diameter_max, diameter_min);
25
26 if (!(*grains = calloc(n[0] * n[1] * n[2], sizeof(struct grain))))
27 err(1, "%s: grains calloc", __func__);
28
29 for (k = 0; k < n[2]; k++)
30 for (j = 0; j < n[1]; j++)
31 for (i = 0; i < n[0]; i++) {
32 N++;
33 ng = idx3(i, j, k, n[0], n[1]);
34 grain_defaults(&(*grains)[ng]);
35 (*grains)[ng].diameter = gsd(diameter_min, diameter_max);
36 (*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
37 (*grains)[ng].pos[1] = j * dx + 0.5 * dx + origo[1];
38 (*grains)[ng].pos[2] = k * dx + 0.5 * dx + origo[2];
39 for (l = 0; l < 3; l++)
40 (*grains)[ng].pos[l] +=
41 random_value_uniform(-0.5 * dx_padding,
42 0.5 * dx_padding);
43 #ifdef VERBOSE
44 printf("added grain %zu\n", ng);
45 #endif
46 }
47
48 return N;
49 }
50
51 size_t
52 triangular_packing(struct grain **grains,
53 size_t n[3],
54 double diameter_min, double diameter_max,
55 double (*gsd)(double min, double max),
56 double padding_factor,
57 double origo[3])
58 {
59 size_t i, j, k, l, ng, N = 0;
60 double dx_padding = diameter_max * padding_factor;
61 double dx = diameter_max + dx_padding;
62 double dy = dx * sin(PI/3.0);
63
64 if (diameter_max < diameter_min)
65 errx(1, "%s: diameter_max (%g) is smaller than diameter_min (%g)",
66 __func__, diameter_max, diameter_min);
67
68 if (!(*grains = calloc(n[0] * n[1] * n[2], sizeof(struct grain))))
69 err(1, "%s: grains calloc", __func__);
70
71 for (k = 0; k < n[2]; k++)
72 for (j = 0; j < n[1]; j++)
73 for (i = 0; i < n[0]; i++) {
74 N++;
75 ng = idx3(i, j, k, n[0], n[1]);
76 grain_defaults(&(*grains)[ng]);
77 (*grains)[ng].diameter = gsd(diameter_min, diameter_max);
78 (*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
79 if (j % 2 == 0)
80 (*grains)[ng].pos[0] += 0.5 * dx;
81 (*grains)[ng].pos[1] = j * dy + 0.5 * dx + origo[1];
82 (*grains)[ng].pos[2] = k * dx + 0.5 * dx + origo[2];
83 for (l = 0; l < 2; l++)
84 (*grains)[ng].pos[l] +=
85 random_value_uniform(-0.5 * dx_padding,
86 0.5 * dx_padding);
87 #ifdef VERBOSE
88 printf("added grain %zu\n", ng);
89 #endif
90 }
91
92 return N;
93 }