| tarrays.c - cngf-pf - continuum model for granular flows with pore-pressure dyn… | |
| git clone git://src.adamsgaard.dk/cngf-pf | |
| Log | |
| Files | |
| Refs | |
| README | |
| 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 … | |
| 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 } |