From 56f643759e666b2f1bbb348a0aa0c65598614384 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sun, 13 Sep 2009 10:10:18 -0400 Subject: [PATCH] Added full_chain_stiffness() to handle stiffness with piston tension. Currently just sets a hardcoded stiffness dx, but that's pretty ugly. A slightly better solution would be a user-controlled option, but that is low on my agenda... --- src/sawsim.nw | 93 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 89 insertions(+), 4 deletions(-) diff --git a/src/sawsim.nw b/src/sawsim.nw index c247d7c..87ed220 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -2316,6 +2316,7 @@ the implementation in [[tension_balance.c]], and the unit testing in <> <> <> +<> <> @ @@ -2431,9 +2432,14 @@ double tension_balance(int num_tension_groups, F = (*tension_handlers[0])(xo, params[0]); - if (stiffness != NULL) + if (stiffness != NULL) { *stiffness = chain_stiffness(&x_xo_data, min_relx); - + if (*stiffness == 0.0) { /* re-calculate based on full chain */ + *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers, + inverse_tension_handlers, params, + last_x, x, min_relx, F); + } + } return F; } @@ -2499,6 +2505,10 @@ double x_of_xo(double xo, void *pdata) } @ +Determine the stiffness of a single tension group by taking a small +step [[dx]] back from the position [[x]] for which we want the +stiffness. The touchy part is determining a good [[dx]] and ensuring +the whole interval is on [[x>=0]]. <>= double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx) { @@ -2524,10 +2534,19 @@ double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double } @ +Attempt to calculate the chain stiffness by adding group stiffnesses +as springs in series. In the case where a group has undefined +stiffness (e.g. piston tension, Section \ref{sec.piston_tension}), +this algorithm breaks down. In that case, [[tension_balance()]] +(\ref{sec.tension_balance}) should use [[full_chain_stiffness()]] +which uses the full chain's $F(x)$ rather than that of the individual +domains'. We attempt the series approach first because, lacking the +numerical inversion inside [[tension_balance()]], it is both faster +and more accurate than the full-chain derivative. <>= double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx) { - double stiffness; + double stiffness, gstiff; int i; /* add group stiffnesses in series */ for (i=0; i < x_xo_data->n_groups; i++) { @@ -2535,7 +2554,10 @@ double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx) fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx); fflush(stderr); #endif - stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx); + gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx); + if (gstiff == 0.0); + return 0.0; + stiffness += 1.0/gstiff; } stiffness = 1.0 / stiffness; return stiffness; @@ -2543,6 +2565,63 @@ double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx) @ +Determine the chain stiffness using only [[tension_balance()]]. This +works around difficulties with tension models that have undefined +stiffness (see discussion for [[chain_stiffness()]] above). +<>= +double full_chain_stiffness(int num_tension_groups, + one_dim_fn_t **tension_handlers, + one_dim_fn_t **inverse_tension_handlers, + void **params, /* array of pointers to tension_handler_data_t */ + double last_x, + double x, + double relx, + double F /* already evaluated F(x) */) +{ + double dx, xl, Fl, stiffness; + + assert(x >= 0); + assert(relx > 0 && relx < 1); + + /* Other option for dx that we ignore because of our poor tension_balance() + * resolution (i.e. bad slopes if you zoom in too close): + * if (last_x != -1.0 && last_x != x) // last_x limited + * dx fabs(x-last_x); + * else + * dx = HUGE_VAL; + * if (1==1) { // constant max-value limited + * dx_p = 1e-12; + * dx = MIN(dx, dx_p); + * } + * if (x != 0 && relx != 0) { // relx limited + * dx_p = x*relx; + * dx = MIN(dx, dx_p); + * } + * TODO, determine which of (min_relx, min_rely, max_steps) is actually + * limiting tension_balance. + */ + dx = 1e-12; /* HACK, how to get length scale? */ + + xl = x-dx; + if (xl >= 0) { + Fl = tension_balance(num_tension_groups, tension_handlers, + inverse_tension_handlers, params, + last_x, xl, NULL); + } else { + xl = x; + Fl = F; + x += dx; + F = tension_balance(num_tension_groups, tension_handlers, + inverse_tension_handlers, params, + last_x, x, NULL); + } + stiffness = (F-Fl)/dx; + assert(stiffness >= 0); + return stiffness; +} + +@ + Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$. Simple bisection, so it's robust and converges fairly quickly. You can set @@ -4361,6 +4440,12 @@ tension at $x=L$ is undefined, but may be any positive value. The model is called the ``piston'' model because it resembles a frictionless piston sliding in a rigid cylinder of length $L$. +Note that because of it's infinte stiffness at $x=L$, fully extended +domains with this tension model will be identical to completely rigid +domains. However, there is a distinction between this model and rigid +domains for $x