tarrays.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tarrays.c (6437B)
---
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <err.h>
5
6 #define DEG2RAD(x) (x*M_PI/180.0)
7
8 void
9 check_magnitude(const char *func_name, int limit, int value)
10 {
11 if (value < limit)
12 errx(1, "%s: input size %d is less than %d\n",
13 func_name, value, limit);
14 }
15
16 /* Translate a i,j,k index in grid with dimensions nx, ny, nz into a
17 * linear index */
18 unsigned int
19 idx3(const unsigned int i,
20 const unsigned int j,
21 const unsigned int k,
22 const unsigned int nx,
23 const unsigned int ny)
24 {
25 return i + nx * j + nx * ny * k;
26 }
27
28 /* Translate a i,j,k index in grid with dimensions nx, ny, nz and a
29 * padding of single ghost nodes into a linear index */
30 unsigned int
31 idx3g(const unsigned int i,
32 const unsigned int j,
33 const unsigned int k,
34 const unsigned int nx,
35 const unsigned int ny)
36 {
37 return i + 1 + (nx + 2) * (j + 1) + (nx + 2) * (ny + 2) * (k + 1);
38 }
39
40 /* Translate a i,j index in grid with dimensions nx, ny into a linear
41 * index */
42 unsigned int
43 idx2(const unsigned int i, const unsigned int j, const unsigned int nx)
44 {
45 return i + nx * j;
46 }
47
48 /* Translate a i,j index in grid with dimensions nx, ny and a padding of
49 * single ghost nodes into a linear index */
50 unsigned int
51 idx2g(const unsigned int i, const unsigned int j, const unsigned int nx)
52 {
53 return i + 1 + (nx + 2) * (j + 1);
54 }
55
56 /* Translate a i index in grid with a padding of single into a linear
57 * index */
58 unsigned int
59 idx1g(const unsigned int i)
60 {
61 return i + 1;
62 }
63
64 /* Return an array of `n` linearly spaced values in the range [lower;
65 * upper] */
66 double *
67 linspace(const double lower, const double upper, const int n)
68 {
69 int i;
70 double *x;
71 double dx;
72
73 check_magnitude(__func__, 1, n);
74 if (!(x = calloc(n, sizeof(double))))
75 err(1, "%s: x calloc", __func__);
76 dx = (upper - lower) / (double) (n - 1);
77 for (i = 0; i < n; ++i)
78 x[i] = lower + dx * i;
79
80 return x;
81 }
82
83 /* Return an array of `n-1` values with the intervals between `x` values */
84 double *
85 spacing(const double *x, const int n)
86 {
87 int i;
88 double *dx;
89
90 check_magnitude(__func__, 2, n);
91 if (!(dx = calloc((n - 1), sizeof(double))))
92 err(1, "%s: dx calloc", __func__);
93 for (i = 0; i < n - 1; ++i)
94 dx[i] = x[i + 1] - x[i];
95
96 return dx;
97 }
98
99 /* Return an array of `n` values with the value 0.0 */
100 double *
101 zeros(const int n)
102 {
103 int i;
104 double *x;
105
106 check_magnitude(__func__, 1, n);
107 if (!(x = calloc(n, sizeof(double))))
108 err(1, "%s: x calloc", __func__);
109 for (i = 0; i < n; ++i)
110 x[i] = 0.0;
111
112 return x;
113 }
114
115 /* Return an array of `n` values with the value 1.0 */
116 double *
117 ones(const int n)
118 {
119 int i;
120 double *x;
121
122 check_magnitude(__func__, 1, n);
123 if (!(x = calloc(n, sizeof(double))))
124 err(1, "%s: x calloc", __func__);
125 for (i = 0; i < n; ++i)
126 x[i] = 1.0;
127
128 return x;
129 }
130
131 /* Return an array of `n` values with a specified value */
132 double *
133 initval(const double value, const int n)
134 {
135 int i;
136 double *x;
137
138 check_magnitude(__func__, 1, n);
139 if (!(x = calloc(n, sizeof(double))))
140 err(1, "%s: x calloc", __func__);
141 for (i = 0; i < n; ++i)
142 x[i] = value;
143
144 return x;
145 }
146
147 /* Return an array of `n` uninitialized values */
148 double *
149 empty(const int n)
150 {
151 double *out;
152
153 check_magnitude(__func__, 1, n);
154
155 if (!(out = calloc(n, sizeof(double))))
156 err(1, "%s: calloc", __func__);
157
158 return out;
159 }
160
161 /* Return largest value in array `a` with size `n` */
162 double
163 max(const double *a, const int n)
164 {
165 int i;
166 double maxval;
167
168 check_magnitude(__func__, 1, n);
169 maxval = -INFINITY;
170 for (i = 0; i < n; ++i)
171 if (a[i] > maxval)
172 maxval = a[i];
173
174 return maxval;
175 }
176
177 /* Return smallest value in array `a` with size `n` */
178 double
179 min(const double *a, const int n)
180 {
181 int i;
182 double minval;
183
184 check_magnitude(__func__, 1, n);
185 minval = +INFINITY;
186 for (i = 0; i < n; ++i)
187 if (a[i] < minval)
188 minval = a[i];
189
190 return minval;
191 }
192
193 void
194 print_array(const double *a, const int n)
195 {
196 int i;
197
198 check_magnitude(__func__, 1, n);
199 for (i = 0; i < n; ++i)
200 printf("%.17g\n", a[i]);
201 }
202
203 void
204 print_arrays(const double *a, const double *b, const int n)
205 {
206 int i;
207
208 check_magnitude(__func__, 1, n);
209 for (i = 0; i < n; ++i)
210 printf("%.17g\t%.17g\n", a[i], b[i]);
211 }
212
213 void
214 print_arrays_2nd_normalized(const double *a, const double *b, const int n)
215 {
216 int i;
217 double max_b;
218
219 check_magnitude(__func__, 1, n);
220 max_b = max(b, n);
221 for (i = 0; i < n; ++i)
222 printf("%.17g\t%.17g\n", a[i], b[i] / max_b);
223 }
224
225 void
226 print_three_arrays(const double *a,
227 const double *b,
228 const double *c,
229 const int n)
230 {
231 int i;
232
233 check_magnitude(__func__, 1, n);
234 for (i = 0; i < n; ++i)
235 printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
236 }
237
238 void
239 fprint_arrays(FILE *fp, const double *a, const double *b, const int n)
240 {
241 int i;
242
243 check_magnitude(__func__, 1, n);
244 for (i = 0; i < n; ++i)
245 fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]);
246 }
247
248 void
249 fprint_three_arrays(FILE *fp,
250 const double *a,
251 const double *b,
252 const double *c,
253 const int n)
254 {
255 int i;
256
257 check_magnitude(__func__, 1, n);
258 for (i = 0; i < n; ++i)
259 fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
260 }
261
262 void
263 copy_values(const double *in, double *out, const int n)
264 {
265 int i;
266
267 check_magnitude(__func__, 1, n);
268 for (i = 0; i < n; ++i)
269 out[i] = in[i];
270 }
271
272 double *
273 copy(const double *in, const int n)
274 {
275 double *out;
276
277 check_magnitude(__func__, 1, n);
278 out = empty(n);
279 copy_values(in, out, n);
280 return out;
281 }
282
283 double *
284 normalize(const double *in, const int n)
285 {
286 int i;
287 double max_val;
288 double *out;
289
290 check_magnitude(__func__, 1, n);
291 out = empty(n);
292 copy_values(in, out, n);
293 max_val = max(out, n);
294
295 if (max_val == 0.0)
296 max_val = 1.0;
297
298 for (i = 0; i < n; ++i)
299 out[i] /= max_val;
300
301 return out;
302 }
303
304 double
305 euclidean_norm(const double *a, const int n)
306 {
307 int i;
308 double out = 0.0;
309
310 for (i = 0; i < n; i++)
311 out += pow(a[i], 2.0);
312
313 return sqrt(out);
314 }
315
316 double
317 euclidean_distance(const double *a, const double *b, const int n)
318 {
319 int i;
320 double out = 0.0;
321
322 for (i = 0; i < n; i++)
323 out += pow(b[i] - a[i], 2.0);
324
325 return sqrt(out);
326 }
327
328 double
329 dot(const double *a, const double *b, const int n)
330 {
331 int i;
332 double out = 0.0;
333
334 for (i = 0; i < n; i++)
335 out += a[i] * b[i];
336
337 return out;
338 }
339
340 double *
341 cross(const double a[3], const double b[3])
342 {
343 double *out = empty(3);
344
345 out[0] = a[1]*b[2] - a[2]*b[1];
346 out[1] = a[2]*b[0] - a[0]*b[2];
347 out[2] = a[0]*b[1] - a[1]*b[0];
348
349 return out;
350 }