tgrain.c - granular - granular dynamics simulation (HTM) git clone git://src.adamsgaard.dk/granular (DIR) Log (DIR) Files (DIR) Refs (DIR) README (DIR) LICENSE --- tgrain.c (11985B) --- 1 #include <stdio.h> 2 #include <stdlib.h> 3 #include <math.h> 4 #include <err.h> 5 #include "granular.h" 6 #include "util.h" 7 #include "grain.h" 8 #include "arrays.h" 9 #include "contact.h" 10 11 #define FLOATPREC 17 12 13 struct grain 14 grain_new(void) 15 { 16 struct grain g; 17 grain_defaults(&g); 18 return g; 19 } 20 21 void 22 grain_defaults(struct grain *g) 23 { 24 int d; 25 26 g->diameter = 1.0; 27 for (d = 0; d < 3; d++) { 28 g->pos[d] 29 = g->vel[d] 30 = g->acc[d] 31 = g->force[d] 32 = g->angpos[d] 33 = g->angvel[d] 34 = g->angacc[d] 35 = g->torque[d] 36 = g->disp[d] 37 = g->forceext[d] 38 = g->contact_stress[d] 39 = 0.0; 40 g->gridpos[d] 41 = g->acc_lock[d] = 0; 42 } 43 g->density = 2.5e3; 44 g->fixed = 0; 45 g->rotating = 1; 46 g->enabled = 1; 47 g->youngs_modulus = 1e9; 48 g->poissons_ratio = 0.2; 49 g->friction_coeff = 0.4; 50 g->damping_n = 0.0; 51 g->damping_t = 0.0; 52 g->tensile_strength = 1e8; 53 g->shear_strength = 1e8; 54 g->fracture_toughness = 1e8; 55 g->ncontacts = 0; 56 for (d = 0; d < MAXCONTACTS; d++) 57 contact_defaults(&g->contacts[d]); 58 g->thermal_energy = 0.0; 59 g->color = 0; 60 } 61 62 static void 63 print_padded_nd_double(FILE *stream, const double *arr) 64 { 65 int d; 66 67 for (d = 0; d < 3; d++) 68 fprintf(stream, "%.*g\t", FLOATPREC, arr[d]); 69 } 70 71 static void 72 print_padded_nd_int(FILE *stream, const int *arr) 73 { 74 int d; 75 76 for (d = 0; d < 3; d++) 77 fprintf(stream, "%d\t", arr[d]); 78 } 79 80 static void 81 print_padded_nd_uint(FILE *stream, const size_t *arr) 82 { 83 int d; 84 85 for (d = 0; d < 3; d++) 86 fprintf(stream, "%zu\t", arr[d]); 87 } 88 89 void 90 grain_print(FILE *stream, const struct grain *g) 91 { 92 fprintf(stream, "%.*g\t", FLOATPREC, g->diameter); 93 print_padded_nd_double(stream, g->pos); 94 print_padded_nd_double(stream, g->vel); 95 print_padded_nd_double(stream, g->acc); 96 print_padded_nd_int(stream, g->acc_lock); 97 print_padded_nd_double(stream, g->force); 98 print_padded_nd_double(stream, g->angpos); 99 print_padded_nd_double(stream, g->angvel); 100 print_padded_nd_double(stream, g->angacc); 101 print_padded_nd_double(stream, g->torque); 102 print_padded_nd_double(stream, g->disp); 103 print_padded_nd_double(stream, g->forceext); 104 fprintf(stream, "%.*g\t", FLOATPREC, g->density); 105 fprintf(stream, "%d\t", g->fixed); 106 fprintf(stream, "%d\t", g->rotating); 107 fprintf(stream, "%d\t", g->enabled); 108 fprintf(stream, "%.*g\t", FLOATPREC, g->youngs_modulus); 109 fprintf(stream, "%.*g\t", FLOATPREC, g->poissons_ratio); 110 fprintf(stream, "%.*g\t", FLOATPREC, g->friction_coeff); 111 fprintf(stream, "%.*g\t", FLOATPREC, g->damping_n); 112 fprintf(stream, "%.*g\t", FLOATPREC, g->damping_t); 113 fprintf(stream, "%.*g\t", FLOATPREC, g->tensile_strength); 114 fprintf(stream, "%.*g\t", FLOATPREC, g->shear_strength); 115 fprintf(stream, "%.*g\t", FLOATPREC, g->fracture_toughness); 116 print_padded_nd_uint(stream, g->gridpos); 117 fprintf(stream, "%zu\t", g->ncontacts); 118 print_padded_nd_double(stream, g->contact_stress); 119 fprintf(stream, "%.*g\t", FLOATPREC, g->thermal_energy); 120 fprintf(stream, "%d\n", g->color); 121 } 122 123 struct grain * 124 grain_read(char *line) 125 { 126 struct grain *g; 127 128 if (!(g = malloc(sizeof(struct grain)))) 129 err(1, "%s: grain malloc", __func__); 130 131 if (sscanf(line, "%lg\t" /* diameter */ 132 "%lg\t%lg\t%lg\t" /* pos */ 133 "%lg\t%lg\t%lg\t" /* vel */ 134 "%lg\t%lg\t%lg\t" /* acc */ 135 "%d\t%d\t%d\t" /* acc_lock */ 136 "%lg\t%lg\t%lg\t" /* force */ 137 "%lg\t%lg\t%lg\t" /* angpos */ 138 "%lg\t%lg\t%lg\t" /* angvel */ 139 "%lg\t%lg\t%lg\t" /* angacc */ 140 "%lg\t%lg\t%lg\t" /* torque */ 141 "%lg\t%lg\t%lg\t" /* disp */ 142 "%lg\t%lg\t%lg\t" /* forceext */ 143 "%lg\t" /* density */ 144 "%d\t" /* fixed */ 145 "%d\t" /* rotating */ 146 "%d\t" /* enabled */ 147 "%lg\t" /* youngs_modulus */ 148 "%lg\t" /* poissons_ratio */ 149 "%lg\t" /* friction_coeff */ 150 "%lg\t" /* damping_n */ 151 "%lg\t" /* damping_t */ 152 "%lg\t" /* tensile_strength */ 153 "%lg\t" /* shear_strength */ 154 "%lg\t" /* fracture_toughness */ 155 "%zu\t%zu\t%zu\t" /* gridpos */ 156 "%zu\t" /* ncontacts */ 157 "%lg\t%lg\t%lg\t" /* contact_stress */ 158 "%lg\t" /* thermal_energy */ 159 "%d\n", /* color */ 160 &g->diameter, 161 &g->pos[0], &g->pos[1], &g->pos[2], 162 &g->vel[0], &g->vel[1], &g->vel[2], 163 &g->acc[0], &g->acc[1], &g->acc[2], 164 &g->acc_lock[0], &g->acc_lock[1], &g->acc_lock[2], 165 &g->force[0], &g->force[1], &g->force[2], 166 &g->angpos[0], &g->angpos[1], &g->angpos[2], 167 &g->angvel[0], &g->angvel[1], &g->angvel[2], 168 &g->angacc[0], &g->angacc[1], &g->angacc[2], 169 &g->torque[0], &g->torque[1], &g->torque[2], 170 &g->disp[0], &g->disp[1], &g->disp[2], 171 &g->forceext[0], &g->forceext[1], &g->forceext[2], 172 &g->density, 173 &g->fixed, 174 &g->rotating, 175 &g->enabled, 176 &g->youngs_modulus, 177 &g->poissons_ratio, 178 &g->friction_coeff, 179 &g->damping_n, 180 &g->damping_t, 181 &g->tensile_strength, 182 &g->shear_strength, 183 &g->fracture_toughness, 184 &g->gridpos[0], &g->gridpos[1], &g->gridpos[2], 185 &g->ncontacts, 186 &g->contact_stress[0], &g->contact_stress[1], &g->contact_stress[2], 187 &g->thermal_energy, 188 &g->color) != 55) 189 errx(1, "%s: could not read line: %s", __func__, line); 190 191 if (grain_check_values(g)) 192 errx(1, "%s: invalid values in grain line: %s", __func__, line); 193 194 return g; 195 } 196 197 int 198 grain_check_values(const struct grain *g) 199 { 200 int d; 201 int status = 0; 202 203 check_float_positive("grain->diameter", g->diameter, &status); 204 205 for (d = 0; d < 3; d++) { 206 check_float("grain->pos", g->pos[d], &status); 207 check_float("grain->vel", g->vel[d], &status); 208 check_float("grain->acc", g->acc[d], &status); 209 check_int_bool("grain->acc_lock", g->acc_lock[d], &status); 210 check_float("grain->force", g->force[d], &status); 211 check_float("grain->angpos", g->angpos[d], &status); 212 check_float("grain->angvel", g->angvel[d], &status); 213 check_float("grain->angacc", g->angacc[d], &status); 214 check_float("grain->torque", g->torque[d], &status); 215 check_float("grain->disp", g->disp[d], &status); 216 check_float("grain->forceext", g->forceext[d], &status); 217 check_float("grain->contact_stress", g->contact_stress[d], &status); 218 } 219 220 check_float_non_negative("grain->density", g->density, &status); 221 check_int_bool("grain->fixed", g->fixed, &status); 222 check_int_bool("grain->rotating", g->rotating, &status); 223 check_int_bool("grain->enabled", g->enabled, &status); 224 225 check_float_non_negative("grain->youngs_modulus", 226 g->youngs_modulus, 227 &status); 228 check_float_non_negative("grain->poissons_ratio", 229 g->poissons_ratio, 230 &status); 231 check_float_non_negative("grain->friction_coeff", 232 g->friction_coeff, 233 &status); 234 check_float_non_negative("grain->damping_n", 235 g->damping_n, 236 &status); 237 check_float_non_negative("grain->damping_t", 238 g->damping_t, 239 &status); 240 check_float_non_negative("grain->tensile_strength", 241 g->tensile_strength, 242 &status); 243 check_float_non_negative("grain->shear_strength", 244 g->shear_strength, 245 &status); 246 check_float_non_negative("grain->fracture_toughness", 247 g->fracture_toughness, 248 &status); 249 250 for (d = 0; d < 3; d++) 251 if (g->gridpos[d] > 1) 252 warn_parameter_value("grain->gridpos is not 0 or 1", 253 (double)g->gridpos[d], &status); 254 255 check_float_non_negative("grain->thermal_energy", g->thermal_energy, &status); 256 257 return status; 258 } 259 260 double 261 grain_volume(const struct grain *g) 262 { 263 return 4.0 / 3.0 * PI * pow(g->diameter / 2.0, 3.0); 264 } 265 266 double 267 grain_mass(const struct grain *g) 268 { 269 return grain_volume(g) * g->density; 270 } 271 272 double 273 grain_moment_of_inertia(const struct grain *g) 274 { 275 return 2.0 / 5.0 * grain_mass(g) * pow(g->diameter / 2.0, 2.0); 276 } 277 278 void 279 grain_zero_kinematics(struct grain *g) 280 { 281 int d; 282 283 for (d = 0; d < 3; d++) { 284 g->vel[d] = 0.0; 285 g->acc[d] = 0.0; 286 g->force[d] = 0.0; 287 g->angvel[d] = 0.0; 288 g->angacc[d] = 0.0; 289 g->torque[d] = 0.0; 290 g->disp[d] = 0.0; 291 } 292 } 293 294 void 295 grain_force_reset(struct grain *g) 296 { 297 int d; 298 299 for (d = 0; d < 3; d++) 300 g->force[d] = 0.0; 301 } 302 303 double 304 grain_translational_kinetic_energy(const struct grain *g) 305 { 306 return 0.5 * grain_mass(g) * euclidean_norm(g->vel, 3.0); 307 } 308 309 double 310 grain_rotational_kinetic_energy(const struct grain *g) 311 { 312 return 0.5 * grain_moment_of_inertia(g) * euclidean_norm(g->angvel, 3.0); 313 } 314 315 double 316 grain_kinetic_energy(const struct grain *g) 317 { 318 return grain_translational_kinetic_energy(g) + 319 grain_rotational_kinetic_energy(g); 320 } 321 322 double 323 grain_thermal_energy(const struct grain *g) 324 { 325 return g->thermal_energy; 326 } 327 328 static void 329 grain_temporal_integration_two_term_taylor(struct grain *g, double dt) 330 { 331 int d; 332 double dx, mass = grain_mass(g); 333 double moment_of_inertia = grain_moment_of_inertia(g); 334 335 for (d = 0; d < 3; d++) { 336 g->acc[d] = (g->force[d] + g->forceext[d]) / mass; 337 338 if (g->rotating) 339 g->angacc[d] = g->torque[d] / moment_of_inertia; 340 } 341 342 if (g->fixed) 343 for (d = 0; d < 3; d++) { 344 g->angacc[d] = 0.0; 345 if (!g->acc_lock[d]) 346 g->acc[d] = 0.0; 347 } 348 349 for (d = 0; d < 3; d++) { 350 dx = g->vel[d] * dt + 0.5 * g->acc[d] * dt * dt; 351 g->pos[d] += dx; 352 g->disp[d] += dx; 353 g->angpos[d] += g->angvel[d] * dt + 0.5 * g->angacc[d] * dt * dt; 354 g->vel[d] += g->acc[d] * dt; 355 g->angvel[d] += g->angacc[d] * dt; 356 } 357 } 358 359 void 360 grain_temporal_integration(struct grain *g, double dt) 361 { 362 grain_temporal_integration_two_term_taylor(g, dt); 363 } 364 365 void 366 grain_register_contact(struct grain *g, size_t i, size_t j, 367 double centerdist[3], double overlap) 368 { 369 int ic, d; 370 371 /* first pass: check if contact is already registered */ 372 for (ic = 0; ic < MAXCONTACTS; ic++) 373 if (g->contacts[ic].active && g->contacts[ic].j == j) { 374 g->contacts[ic].overlap = overlap; 375 for (d = 0; d < 3; d++) 376 g->contacts[ic].centerdist[d] = centerdist[d]; 377 return; 378 } 379 380 /* second pass: register as new contact if overlapping */ 381 if (overlap > 0.0) 382 for (ic = 0; ic <= MAXCONTACTS; ic++) { 383 if (ic == MAXCONTACTS) 384 err(1, "%s: contact %zu exceeds MAXCONTACTS (%d)", 385 __func__, i, MAXCONTACTS); 386 if (!g->contacts[ic].active) { 387 contact_new(&g->contacts[ic], i, j); 388 g->ncontacts += 1; 389 g->contacts[ic].overlap = overlap; 390 for (d = 0; d < 3; d++) { 391 g->contacts[ic].centerdist[d] = centerdist[d]; 392 g->contacts[ic].tandisp[d] = 0.0; 393 } 394 return; 395 } 396 } 397 } 398 399 double 400 grain_stiffness_normal(const struct grain *g) 401 { 402 return 2.0 / 3.0 * g->youngs_modulus / (1.0 - pow(g->poissons_ratio, 2.0)); 403 } 404 405 double 406 grain_stiffness_tangential(const struct grain *g) 407 { 408 return 2.0 * g->youngs_modulus 409 / ((1.0 + g->poissons_ratio) * (2.0 - g->poissons_ratio)); 410 } 411 412 double 413 grain_collision_time(const struct grain *g) 414 { 415 double m; 416 417 m = grain_mass(g); 418 419 return fmin(PI * pow(2.0 * grain_stiffness_normal(g) / m 420 - pow(g->damping_n, 2.0) / 4.0, -0.5), 421 PI * pow(2.0 * grain_stiffness_tangential(g) / m 422 - pow(g->damping_t, 2.0) / 4.0, -0.5)); 423 }