tslider.c - slidergrid - grid of elastic sliders on a frictional surface
HTML git clone git://src.adamsgaard.dk/slidergrid
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tslider.c (23400B)
---
1 #include <stdio.h>
2 #include <math.h>
3 #include "typedefs.h"
4 #include "slider.h"
5 #include "vector_math.h"
6 #include "debug.h"
7
8 void print_slider_position(slider s)
9 {
10 printf("%f\t%f\t%f\n", s.pos.x, s.pos.y, s.pos.z);
11 }
12
13 void initialize_slider_values(slider* s)
14 {
15 s->pos_future = zeroes_float3();
16 s->vel = zeroes_float3();
17 s->acc = zeroes_float3();
18 s->angpos_future = zeroes_float3();
19 s->angpos = zeroes_float3();
20 s->angvel = zeroes_float3();
21 s->angacc = zeroes_float3();
22
23 s->force = zeroes_float3();
24 s->torque = zeroes_float3();
25
26 s->mass = 0.0;
27 s->moment_of_inertia = 0.0;
28
29 s->bond_parallel_kv_stiffness = 0.0;
30 s->bond_parallel_kv_viscosity = 0.0;
31
32 s->bond_shear_kv_stiffness = 0.0;
33 s->bond_shear_kv_viscosity = 0.0;
34
35 s->bond_twist_kv_stiffness = 0.0;
36 s->bond_twist_kv_viscosity = 0.0;
37
38 s->bond_bend_kv_stiffness = 0.0;
39 s->bond_bend_kv_viscosity = 0.0;
40
41 s->damping_viscosity_linear = 0.0;
42 s->damping_viscosity_angular = 0.0;
43 s->damping_coefficient = 0.0;
44
45 // define all entries in neighbor list as empty
46 int i;
47 for (i=0; i<MAX_NEIGHBORS; i++) {
48 s->neighbors[i] = -1;
49 s->neighbor_distance[i] = zeroes_float3();
50 s->neighbor_relative_distance_displacement[i] = 0.0;
51 s->neighbor_relative_distance_velocity[i] = zeroes_float3();
52 s->neighbor_relative_tangential_displacement[i] = zeroes_float3();
53 s->neighbor_relative_tangential_velocity[i] = zeroes_float3();
54 s->neighbor_relative_twist[i] = 0.0;
55 s->neighbor_relative_twist_velocity[i] = 0.0;
56 s->neighbor_relative_bend[i] = zeroes_float3();
57 s->neighbor_relative_bend_velocity[i] = zeroes_float3();
58 }
59 }
60
61
62 // find change in inter-slider distance from old position and projected new
63 // position, based on current kinematics and explicit temporal integration using
64 // TY2 expansion
65 void project_slider_position(slider* s, Float dt, long int iteration)
66 {
67 s->pos_future = make_float3(
68 s->pos.x + s->vel.x*dt + 0.5*s->acc.x*dt*dt,
69 s->pos.y + s->vel.y*dt + 0.5*s->acc.y*dt*dt,
70 s->pos.z + s->vel.z*dt + 0.5*s->acc.z*dt*dt);
71
72 s->angpos_future = make_float3(
73 s->angpos.x + s->angvel.x*dt + 0.5*s->angacc.x*dt*dt,
74 s->angpos.y + s->angvel.y*dt + 0.5*s->angacc.y*dt*dt,
75 s->angpos.z + s->angvel.z*dt + 0.5*s->angacc.z*dt*dt);
76 }
77
78 // reset sum of forces and torques on slider
79 void zero_slider_force_torque(slider* s)
80 {
81 s->force = zeroes_float3();
82 s->torque = zeroes_float3();
83 }
84
85 /* Explicit temporal integration scheme based on three-term Taylor expansion.
86 * Truncation error O(dt^4) for positions, O(dt^3) for velocities. Acceleration
87 * change is approximated by backwards differences. */
88 void update_kinematics(slider* s, Float dt, long int iteration)
89 {
90 s->acc = divide_float3_scalar(s->force, s->mass);
91 s->angacc = divide_float3_scalar(s->torque, s->moment_of_inertia);
92
93 Float3 acc0, angacc0;
94 if (iteration == 0) {
95 acc0 = zeroes_float3();
96 angacc0 = zeroes_float3();
97 } else {
98 acc0 = s->acc;
99 angacc0 = s->angacc;
100 }
101
102 const Float3 dacc_dt = make_float3(
103 (s->acc.x - acc0.x)/dt,
104 (s->acc.y - acc0.y)/dt,
105 (s->acc.z - acc0.z)/dt);
106
107 const Float3 dangacc_dt = make_float3(
108 (s->angacc.x - angacc0.x)/dt,
109 (s->angacc.y - angacc0.y)/dt,
110 (s->angacc.z - angacc0.z)/dt);
111
112 const Float3 dpos_dt = make_float3(
113 s->vel.x*dt + 0.5*s->acc.x*dt*dt + 1./6.*dacc_dt.x*dt*dt*dt,
114 s->vel.y*dt + 0.5*s->acc.y*dt*dt + 1./6.*dacc_dt.y*dt*dt*dt,
115 s->vel.z*dt + 0.5*s->acc.z*dt*dt + 1./6.*dacc_dt.z*dt*dt*dt);
116
117 const Float3 dangpos_dt = make_float3(
118 s->angvel.x*dt + 0.5*s->angacc.x*dt*dt
119 + 1./6.*dangacc_dt.x*dt*dt*dt,
120 s->angvel.y*dt + 0.5*s->angacc.y*dt*dt
121 + 1./6.*dangacc_dt.y*dt*dt*dt,
122 s->angvel.z*dt + 0.5*s->angacc.z*dt*dt
123 + 1./6.*dangacc_dt.z*dt*dt*dt);
124
125 const Float3 dvel_dt = make_float3(
126 s->acc.x*dt + 0.5*dacc_dt.x*dt*dt,
127 s->acc.y*dt + 0.5*dacc_dt.y*dt*dt,
128 s->acc.z*dt + 0.5*dacc_dt.z*dt*dt);
129
130 const Float3 dangvel_dt = make_float3(
131 s->angacc.x*dt + 0.5*dangacc_dt.x*dt*dt,
132 s->angacc.y*dt + 0.5*dangacc_dt.y*dt*dt,
133 s->angacc.z*dt + 0.5*dangacc_dt.z*dt*dt);
134
135 s->pos = add_float3(s->pos, dpos_dt);
136 s->angpos = add_float3(s->angpos, dangpos_dt);
137 s->vel = add_float3(s->vel, dvel_dt);
138 s->angvel = add_float3(s->angvel, dangvel_dt);
139 }
140
141
142 // determine time step increment bond-parallel deformation (tension or
143 // compression)
144 void bond_tensile_deformation(slider* s1, const slider s2,
145 const int idx_neighbor, const int iteration)
146 {
147 // vector pointing from neighbor (s2) position to this slider position (s1)
148 const Float3 dist = subtract_float3(s1->pos, s2.pos);
149
150 // length of inter-slider vector
151 const Float dist_norm = norm_float3(dist);
152
153 // unit vector pointing from neighbor (s2) to slider s1
154 const Float3 n = divide_float3_scalar(dist, dist_norm);
155
156 #ifdef DEBUG_BOND_DEFORMATION
157 printf("\tdist = %f %f %f\n\tdist_norm = %f\n\tn = %f %f %f\n",
158 dist.x,
159 dist.y,
160 dist.z,
161 dist_norm,
162 n.x,
163 n.y,
164 n.z);
165 #endif
166
167 // relative contact interface velocity w/o rolling
168 const Float3 vel_linear = subtract_float3(s1->vel, s2.vel);
169
170 // Normal component of the relative contact interface velocity
171 const Float vel_n = -1.0*dot_float3(vel_linear, n);
172
173 // read previous inter-slider distance vector
174 const Float3 dist0 = s1->neighbor_distance[idx_neighbor];
175
176 // determine projected future displacement
177 const Float3 dist_future = subtract_float3(s1->pos_future, s2.pos_future);
178
179 // increment in inter-slider distance with central differences, divide by
180 // two to get displacement over 1 time step
181 const Float3 ddist = divide_float3_scalar(
182 subtract_float3(dist_future, dist0), 2.0);
183
184 // save current inter-slider distance
185 s1->neighbor_distance[idx_neighbor] = dist;
186
187 // total relative displacement in inter-slider distance
188 if (iteration == 0)
189 s1->neighbor_relative_distance_displacement[idx_neighbor] =
190 dot_float3(n, ddist);
191 else
192 s1->neighbor_relative_distance_displacement[idx_neighbor] +=
193 dot_float3(n, ddist);
194
195 // total relative displacement in inter-slider distance
196 s1->neighbor_relative_distance_velocity[idx_neighbor] =
197 multiply_float3_scalar(n, vel_n);
198 }
199
200
201 // determine time step increment bond-shear deformation
202 void bond_shear_deformation(slider* s1, const slider s2,
203 const int idx_neighbor, const int iteration, const Float dt)
204 {
205 // vector pointing from neighbor (s2) position to this slider position (s1)
206 const Float3 dist = subtract_float3(s1->pos, s2.pos);
207
208 // length of inter-slider vector
209 const Float dist_norm = norm_float3(dist);
210
211 // unit vector pointing from neighbor (s2) to slider s1
212 const Float3 n = divide_float3_scalar(dist, dist_norm);
213
214 // relative contact interface velocity w/o rolling
215 const Float3 vel_linear = subtract_float3(s1->vel, s2.vel);
216
217 // relative contact interface velocity with rolling
218 const Float3 vel =
219 add_float3(
220 vel_linear,
221 add_float3(
222 cross_float3(divide_float3_scalar(dist, 2.0), s1->angvel),
223 cross_float3(divide_float3_scalar(dist, 2.0), s2.angvel)
224 ));
225
226 // Tangential component of the relative contact interface velocity
227 // Hinrichsen and Wolf 2004, eq. 13.9
228 const Float3 vel_t =
229 subtract_float3(vel, multiply_float3_scalar(n, dot_float3(n, vel)));
230
231 // Read previous (uncorrected) tangential displacement vector
232 Float3 tangential_displacement0_uncor;
233 if (iteration == 0)
234 tangential_displacement0_uncor = zeroes_float3();
235 else
236 tangential_displacement0_uncor =
237 s1->neighbor_relative_tangential_displacement[idx_neighbor];
238
239 // Correct previous tangential displacement vector if the contact plane is
240 // reoriented
241 const Float3 tangential_displacement0 =
242 subtract_float3(tangential_displacement0_uncor,
243 multiply_float3(n,
244 multiply_float3(n, tangential_displacement0_uncor)));
245
246 // project future tangential displacement
247 const Float3 dtangential_displacement = multiply_float3_scalar(vel_t, dt);
248 const Float3 tangential_displacement_future =
249 add_float3(tangential_displacement0, dtangential_displacement);
250
251 // get current tangential displacement from central differences
252 const Float3 tangential_displacement = divide_float3_scalar(
253 subtract_float3(tangential_displacement_future,
254 tangential_displacement0), 2.0);
255
256 // total relative displacement in inter-slider distance
257 if (iteration == 0)
258 s1->neighbor_relative_tangential_displacement[idx_neighbor] =
259 tangential_displacement;
260 else
261 s1->neighbor_relative_tangential_displacement[idx_neighbor] =
262 add_float3(
263 s1->neighbor_relative_tangential_displacement[idx_neighbor],
264 dtangential_displacement);
265
266 // total relative displacement in inter-slider distance
267 s1->neighbor_relative_tangential_velocity[idx_neighbor] = vel_t;
268
269 #ifdef DEBUG_BOND_DEFORMATION
270 printf("\t------\n"
271 "\ttan_disp0_u = %f %f %f\n"
272 "\ttan_disp0 = %f %f %f\n"
273 "\ttan_disp = %f %f %f\n"
274 "\tdtan_disp = %f %f %f\n"
275 "\tvel_t = %f %f %f\n"
276 "\t------\n"
277 ,
278 tangential_displacement0_uncor.x,
279 tangential_displacement0_uncor.y,
280 tangential_displacement0_uncor.z,
281 tangential_displacement0.x,
282 tangential_displacement0.y,
283 tangential_displacement0.z,
284 tangential_displacement.x,
285 tangential_displacement.y,
286 tangential_displacement.z,
287 dtangential_displacement.x,
288 dtangential_displacement.y,
289 dtangential_displacement.z,
290 vel_t.x,
291 vel_t.y,
292 vel_t.z);
293 #endif
294 }
295
296 // determine time step increment bond twist deformation
297 void bond_twist_deformation(slider* s1, const slider s2,
298 const int idx_neighbor, const int iteration, const Float dt)
299 {
300 // vector pointing from neighbor (s2) position to this slider position (s1)
301 const Float3 dist = subtract_float3(s1->pos, s2.pos);
302
303 // length of inter-slider vector
304 const Float dist_norm = norm_float3(dist);
305
306 // unit vector pointing from neighbor (s2) to slider s1
307 const Float3 n = divide_float3_scalar(dist, dist_norm);
308
309 // Get the relative angular velocity
310 const Float3 angvel_rel = subtract_float3(s2.angvel, s1->angvel);
311
312 // Isolate the relative angular velocity component parallel to the bond
313 const Float twist_vel = dot_float3(angvel_rel, n);
314
315 // Twist displacement happening in one time step
316 const Float dtwist = twist_vel*dt;
317
318 // total relative displacement in inter-slider distance
319 if (iteration == 0)
320 s1->neighbor_relative_twist[idx_neighbor] = dtwist;
321 else
322 s1->neighbor_relative_twist[idx_neighbor] += dtwist;
323
324 // total relative displacement in inter-slider distance
325 s1->neighbor_relative_twist_velocity[idx_neighbor] = twist_vel;
326
327 #ifdef DEBUG_BOND_DEFORMATION
328 printf("\t------\n"
329 "\ttwist = %f %f %f\n"
330 "\tdtwist = %f %f %f\n"
331 "\ttwist_vel = %f %f %f\n"
332 "\t------\n"
333 ,
334 twist.x,
335 twist.y,
336 twist.z,
337 dtwist.x,
338 dtwist.y,
339 dtwist.z,
340 twist_vel.x,
341 twist_vel.y,
342 twist_vel.z);
343 #endif
344 }
345
346
347
348 // Find the bond deformation
349 void bond_deformation(slider* s1, const slider s2,
350 const int idx_neighbor, const int iteration, const Float dt)
351 {
352 bond_tensile_deformation(s1, s2, idx_neighbor, iteration);
353 bond_shear_deformation(s1, s2, idx_neighbor, iteration, dt);
354 bond_twist_deformation(s1, s2, idx_neighbor, iteration, dt);
355 }
356
357
358 void bond_parallel_kelvin_voigt(slider* s1, const slider s2,
359 const int idx_neighbor)
360 {
361 // unit vector pointing from neighbor (s2) to slider s1
362 const Float3 n = divide_float3_scalar(s1->neighbor_distance[idx_neighbor],
363 norm_float3(s1->neighbor_distance[idx_neighbor]));
364
365 // determine the bond-parallel KV stiffness from the harmonic mean.
366 // differs from Potyondy & Stack 2004, eq. 6.
367 const Float bond_parallel_kv_stiffness =
368 2.*s1->bond_parallel_kv_stiffness*s2.bond_parallel_kv_stiffness/
369 (s1->bond_parallel_kv_stiffness + s2.bond_parallel_kv_stiffness
370 + 1.0e-40);
371
372 // determine the bond-parallel KV viscosity from the harmonic mean.
373 const Float bond_parallel_kv_viscosity =
374 2.*s1->bond_parallel_kv_viscosity*s2.bond_parallel_kv_viscosity/
375 (s1->bond_parallel_kv_viscosity + s2.bond_parallel_kv_viscosity
376 + 1.0e-40);
377
378 // bond-parallel Kelvin-Voigt elasticity
379 const Float3 f_n_elastic =
380 multiply_scalar_float3(bond_parallel_kv_stiffness *
381 s1->neighbor_relative_distance_displacement[idx_neighbor],
382 n);
383
384 // bond-parallel Kelvin-Voigt viscosity
385 const Float3 f_n_viscous =
386 multiply_scalar_float3(bond_parallel_kv_viscosity,
387 s1->neighbor_relative_distance_velocity[idx_neighbor]);
388
389 // bond-parallel Kelvin-Voigt force, counteracts tension and compression
390 const Float3 f_n = multiply_scalar_float3(-1.0,
391 add_float3(f_n_elastic, f_n_viscous));
392
393 // add bond-parallel Kelvin-Voigt force to sum of forces on slider
394 s1->force = add_float3(s1->force, f_n);
395
396 #ifdef DEBUG_SLIDER_FORCE_COMPONENTS
397 printf(" slider_interaction KV-parallel:\n"
398 "\tn = %f %f %f\n\tf_n = %f %f %f\n\t"
399 "f_n_elastic = %f %f %f\n\tf_n_viscous = %f %f %f\n",
400 n.x,
401 n.y,
402 n.z,
403 f_n.x,
404 f_n.y,
405 f_n.z,
406 f_n_elastic.x,
407 f_n_elastic.y,
408 f_n_elastic.z,
409 f_n_viscous.x,
410 f_n_viscous.y,
411 f_n_viscous.z);
412 #endif
413 }
414
415 void bond_shear_kelvin_voigt(slider* s1, const slider s2,
416 const int idx_neighbor)
417 {
418 // determine the bond-shear KV stiffness from the harmonic mean.
419 const Float bond_shear_kv_stiffness =
420 2.*s1->bond_shear_kv_stiffness*s2.bond_shear_kv_stiffness/
421 (s1->bond_shear_kv_stiffness + s2.bond_shear_kv_stiffness
422 + 1.0e-40);
423
424 // determine the bond-shear KV viscosity from the harmonic mean.
425 const Float bond_shear_kv_viscosity =
426 2.*s1->bond_shear_kv_viscosity*s2.bond_shear_kv_viscosity/
427 (s1->bond_shear_kv_viscosity + s2.bond_shear_kv_viscosity
428 + 1.0e-40);
429
430 // bond-normal Kelvin-Voigt elasticity
431 const Float3 f_t_elastic =
432 multiply_scalar_float3(bond_shear_kv_stiffness,
433 s1->neighbor_relative_tangential_displacement[idx_neighbor]);
434
435 // bond-normal Kelvin-Voigt viscosity
436 const Float3 f_t_viscous =
437 multiply_scalar_float3(bond_shear_kv_viscosity,
438 s1->neighbor_relative_tangential_velocity[idx_neighbor]);
439
440 // bond-normal Kelvin-Voigt force, counteracts shear
441 const Float3 f_t = multiply_scalar_float3( -1.0,
442 add_float3(f_t_elastic, f_t_viscous));
443
444 // add bond-shear Kelvin-Voigt force to sum of forces on slider
445 s1->force = add_float3(s1->force, f_t);
446
447 // determine torque on slider from shear on this bond
448 const Float3 torque = multiply_scalar_float3( -1.0,
449 cross_float3(
450 multiply_float3_scalar(
451 s1->neighbor_distance[idx_neighbor], 0.5),
452 f_t));
453
454 // add bond-shear Kelvin-Voigt force to sum of torques on slider
455 s1->torque = add_float3(s1->torque, torque);
456
457 #ifdef DEBUG_SLIDER_FORCE_COMPONENTS
458 printf(" slider_interaction KV-shear:\n"
459 "\tf_t = %f %f %f\n"
460 "\ttorque = %f %f %f\n"
461 "\tf_t_elastic = %f %f %f\n"
462 "\tf_t_viscous = %f %f %f\n"
463 "\ttan_disp = %f %f %f\n"
464 "\ttan_vel = %f %f %f\n",
465 f_t.x,
466 f_t.y,
467 f_t.z,
468 torque.x,
469 torque.y,
470 torque.z,
471 f_t_elastic.x,
472 f_t_elastic.y,
473 f_t_elastic.z,
474 f_t_viscous.x,
475 f_t_viscous.y,
476 f_t_viscous.z,
477 s1->neighbor_relative_tangential_displacement[idx_neighbor].x,
478 s1->neighbor_relative_tangential_displacement[idx_neighbor].y,
479 s1->neighbor_relative_tangential_displacement[idx_neighbor].z,
480 s1->neighbor_relative_tangential_velocity[idx_neighbor].x,
481 s1->neighbor_relative_tangential_velocity[idx_neighbor].y,
482 s1->neighbor_relative_tangential_velocity[idx_neighbor].z);
483 #endif
484 }
485
486 void bond_twist_kelvin_voigt(slider* s1, const slider s2,
487 const int idx_neighbor)
488 {
489 // determine the bond-twist KV stiffness from the harmonic mean.
490 const Float bond_twist_kv_stiffness =
491 2.*s1->bond_twist_kv_stiffness*s2.bond_twist_kv_stiffness/
492 (s1->bond_twist_kv_stiffness + s2.bond_twist_kv_stiffness
493 + 1.0e-40);
494
495 // determine the bond-twist KV viscosity from the harmonic mean.
496 const Float bond_twist_kv_viscosity =
497 2.*s1->bond_twist_kv_viscosity*s2.bond_twist_kv_viscosity/
498 (s1->bond_twist_kv_viscosity + s2.bond_twist_kv_viscosity
499 + 1.0e-40);
500
501 // bond-parallel Kelvin-Voigt elasticity
502 const Float m_n_elastic = bond_twist_kv_stiffness
503 *s1->neighbor_relative_twist[idx_neighbor];
504
505 // bond-parallel Kelvin-Voigt viscosity
506 const Float m_n_viscous = bond_twist_kv_viscosity
507 *s1->neighbor_relative_twist_velocity[idx_neighbor];
508
509 // bond-parallel Kelvin-Voigt moment, counteracts twisting
510 const Float m_n = -1.0*(m_n_elastic + m_n_viscous);
511
512 // determine torque on slider from shear on this bond
513 const Float3 torque = multiply_scalar_float3( -1.0,
514 cross_float3(
515 multiply_float3_scalar(
516 s1->neighbor_distance[idx_neighbor], 0.5),
517 f_t));*/
518
519 // add bond-shear Kelvin-Voigt force to sum of torques on slider
520 s1->torque = add_float3(s1->torque, torque);
521
522 #ifdef DEBUG_SLIDER_FORCE_COMPONENTS
523 printf(" slider_interaction KV-twist:\n"
524 "\tf_t = %f %f %f\n"
525 "\ttorque = %f %f %f\n"
526 "\tf_t_elastic = %f %f %f\n"
527 "\tf_t_viscous = %f %f %f\n"
528 "\ttan_disp = %f %f %f\n"
529 "\ttan_vel = %f %f %f\n",
530 f_t.x,
531 f_t.y,
532 f_t.z,
533 torque.x,
534 torque.y,
535 torque.z,
536 f_t_elastic.x,
537 f_t_elastic.y,
538 f_t_elastic.z,
539 f_t_viscous.x,
540 f_t_viscous.y,
541 f_t_viscous.z,
542 s1->neighbor_relative_tangential_displacement[idx_neighbor].x,
543 s1->neighbor_relative_tangential_displacement[idx_neighbor].y,
544 s1->neighbor_relative_tangential_displacement[idx_neighbor].z,
545 s1->neighbor_relative_tangential_velocity[idx_neighbor].x,
546 s1->neighbor_relative_tangential_velocity[idx_neighbor].y,
547 s1->neighbor_relative_tangential_velocity[idx_neighbor].z);
548 #endif
549 }
550
551 // Resolve bond mechanics between a slider and one of its neighbors based on the
552 // incremental linear-elastic model by Potyondy and Cundall 2004
553 void slider_interaction(slider* s1, const slider s2, const int idx_neighbor)
554 {
555 bond_parallel_kelvin_voigt(s1, s2, idx_neighbor);
556 bond_shear_kelvin_voigt(s1, s2, idx_neighbor);
557 bond_twist_kelvin_voigt(s1, s2, idx_neighbor);
558 }
559
560
561 // Resolve interaction between a slider and all of its neighbors
562 void slider_neighbor_interaction(
563 slider* s,
564 const slider* sliders,
565 const int N,
566 const int iteration,
567 const Float dt)
568 {
569 int idx_neighbor;
570 for (idx_neighbor=0; idx_neighbor<MAX_NEIGHBORS; idx_neighbor++) {
571
572 if (s->neighbors[idx_neighbor] != -1) {
573
574 #ifdef DEBUG_BOND_DEFORMATION
575 printf("%p = %d:\n", s, idx_neighbor);
576 #endif
577 bond_deformation(
578 s, sliders[s->neighbors[idx_neighbor]],
579 idx_neighbor, iteration, dt);
580
581 #ifdef DEBUG_BOND_DEFORMATION
582 printf("\trel_dist_disp = %f \n"
583 "\trel_dist_vel = %f %f %f\n"
584 "\trel_tan_disp = %f %f %f\n"
585 "\trel_tan_vel = %f %f %f\n",
586 s->neighbor_relative_distance_displacement[idx_neighbor],
587 s->neighbor_relative_distance_velocity[idx_neighbor].x,
588 s->neighbor_relative_distance_velocity[idx_neighbor].y,
589 s->neighbor_relative_distance_velocity[idx_neighbor].z,
590 s->neighbor_relative_tangential_displacement[idx_neighbor].x,
591 s->neighbor_relative_tangential_displacement[idx_neighbor].y,
592 s->neighbor_relative_tangential_displacement[idx_neighbor].z,
593 s->neighbor_relative_tangential_velocity[idx_neighbor].x,
594 s->neighbor_relative_tangential_velocity[idx_neighbor].y,
595 s->neighbor_relative_tangential_velocity[idx_neighbor].z
596 );
597 #endif
598
599 slider_interaction(
600 s, sliders[s->neighbors[idx_neighbor]], idx_neighbor);
601
602 /*printf("s.force = %f %f %f\n",
603 s->force.x, s->force.y, s->force.z);*/
604 }
605 }
606 }
607
608 // add frequency-independent viscosity to sum of forces and torques to damp
609 // reflected waves from the edge of the lattice
610 // Mora and Place (1994), Place et al (2002)
611 void slider_viscous_damping(slider* s)
612 {
613 s->force.x += - s->damping_viscosity_linear*s->vel.x;
614 s->force.y += - s->damping_viscosity_linear*s->vel.y;
615 s->force.z += - s->damping_viscosity_linear*s->vel.z;
616
617 s->torque.x += - s->damping_viscosity_angular*s->vel.x;
618 s->torque.y += - s->damping_viscosity_angular*s->vel.y;
619 s->torque.z += - s->damping_viscosity_angular*s->vel.z;
620 }
621
622 // add local non-viscous damping, Potyondy and Cundall (2004)
623 void slider_nonviscous_damping(slider* s)
624 {
625 s->force.x +=
626 - s->damping_coefficient*fabs(s->force.x)*copysign(1.0, s->vel.x);
627 s->force.y +=
628 - s->damping_coefficient*fabs(s->force.y)*copysign(1.0, s->vel.y);
629 s->force.z +=
630 - s->damping_coefficient*fabs(s->force.z)*copysign(1.0, s->vel.z);
631
632 s->torque.x +=
633 - s->damping_coefficient*fabs(s->torque.x)*copysign(1.0, s->angvel.x);
634 s->torque.y +=
635 - s->damping_coefficient*fabs(s->torque.y)*copysign(1.0, s->angvel.y);
636 s->torque.z +=
637 - s->damping_coefficient*fabs(s->torque.z)*copysign(1.0, s->angvel.z);
638
639 }