Mass-Spring-Damper Model
Hi everybody.
I'm trying to figure out a way to implement a physical mass-spring-damper model in Max. Does anybody have an idea how I could do that?
I want to connect 5 faders with imaginary springs with each other. If I move one it should affect the others.
I figured out that there were a bunch of pmpd objects, that did the job, and where you could even define the damping and stiffness of the springs. Unfortunately these object don't seem to work anymore.
Does anybody have an idea?
With kind regards,
Josef
You could take a look here:
http://www.mi-creative.eu/index.html
http://www.mi-creative.eu/post_miTuto_gen.html
mati
hi Josef,
jit.phys is probably your best friend here.
But depending how you want it, you take a look at the "ease" package too.
Otherwise, if you like it to work at audio rate, the best option would be to use a resonant filter.
HTH.
i guess he needs custom (patches) solutions for the WDF thing. ;)
i took a crack at implementing something for this with jit.phys objects. the tricky parts are knowing what attributes of jit.phys.6dof to adjust for spring type constraints, as well as what attributes to 0 out all other constraints aspects, and what settings to tweak to get the right movement.
The attached patch demonstrates one way to do this using a jit.phys.barslide constraint for each slider, and jit.phys.6dof contraints between adjacent bodies for the spring. I found that damping of jit.phys.body and springstiff_lin of jit.phys.6dof are good things to tweak, and I also added a little hack using universal and jit.proxy to take the output of jit.phys.constraint to set the selected body mass to a high value, and all other bodies mass to low whenever a slider is adjusted with the mouse.
we really should come up with a wrapper around 6dof to expose only spring constraints, since this is useful feature and hard to discover without knowing the details of 6dof.
patch attached.
Wow Rob, this is exactly what I had in mind. Great work!
Do you have a suggestion on how I could extract the values, so I can apply them to actual audio faders?
You can ask jit.phys.body's position with the message "getposition", and it will be dumped on the right-most outlet.
@Rob : nice patch! + I agree on the spring constraint being a little bit drowned in all the possibilities offered by 6dof and that a wrapper for that would be nice!
Thanks a lot.
I extended the patch to 16 sliders, and it works exactly how I wanted it.
Now I'm trying to figure out how to 'trigger' a certain fader to go up and down (like I do in the example by hand). Is there a way to do that without including the mouse?
very cool! you probably want to toggle kinematic mode on on the jit.phys.body that you want to manually position, and then toggle it back off when done. this will allow you to move the body via position attribute, and all connected bodies will respond to the simulation forces triggered by the kinematic body.
Thank you Rob, that works. The only thing is that somehow the position of the body rotates 3-dimensionally slowly, as soon as I toggle kinematic mode. In the long term the position gets really funny and creates quite some noise. Any suggestion how I can avoid that?
It looks like this after a while...
yep, sometimes the squares hang on mouse down, i think it looks funny.
otherwise it does exactly the same as the good old pmpd.mass.
which i´ve always would have loved to have in patch form.
Yes a patch of the good old pmpd would be great.
I wander how the mouse down command got in there at the first place? And would you know a way to get rid of it? The whole thing gets really messy after a while like this
https://sources.debian.org/data/main/p/pd-pmpd/0.9-4/manual/pmpd.pdf
page 5.
but still not simple enough for me.
[code]
#include "m_pd.h"
#include "math.h"
#define max(a,b) ( ((a) > (b)) ? (a) : (b) )
#define min(a,b) ( ((a) < (b)) ? (a) : (b) )
static t_class *mass_class;
typedef struct _mass {
t_object x_obj;
t_float pos_old_1, pos_old_2, Xinit;
t_float force, mass, dX;
t_float minX, maxX;
t_outlet *position_new, *vitesse_out, *force_out;
t_symbol *x_sym; // receive
unsigned int x_state; // random
t_float x_f; // random
} t_mass;
static int makeseed(void)
{
static unsigned int random_nextseed = 1489853723;
random_nextseed = random_nextseed * 435898247 + 938284287;
return (random_nextseed & 0x7fffffff);
}
static float random_bang(t_mass *x)
{
int nval;
int range = 2000000;
float rnd;
unsigned int randval = x->x_state;
x->x_state = randval = randval * 472940017 + 832416023;
nval = ((double)range) * ((double)randval)
* (1./4294967296.);
if (nval >= range) nval = range-1;
rnd=nval;
rnd-=1000000;
rnd=rnd/1000000.; //pour mettre entre -1 et 1;
return (rnd);
}
void mass_minX(t_mass *x, t_floatarg f1)
{
x->minX = f1;
}
void mass_maxX(t_mass *x, t_floatarg f1)
{
x->maxX = f1;
}
void mass_float(t_mass *x, t_floatarg f1)
{
x->force += f1;
}
void mass_bang(t_mass *x)
{
t_float pos_new;
if (x->mass > 0)
pos_new = x->force/x->mass + 2*x->pos_old_1 - x->pos_old_2;
else pos_new = x->pos_old_1;
pos_new = max(min(x->maxX, pos_new), x->minX);
pos_new += x->dX;
x->pos_old_1 += x->dX; // pour ne pas avoir d'inertie suplementaire du a ce deplacement
outlet_float(x->vitesse_out, x->pos_old_1 - x->pos_old_2);
outlet_float(x->force_out, x->force);
outlet_float(x->position_new, pos_new);
x->pos_old_2 = x->pos_old_1;
x->pos_old_1 = pos_new;
// x->force = 0;
x->force = random_bang(x)*1e-25; // avoiding denormal problem by adding low amplitude noise
x->dX = 0;
}
void mass_reset(t_mass *x)
{
x->pos_old_2 = x->Xinit;
x->pos_old_1 = x->Xinit;
x->force=0;
outlet_float(x->position_new, x->Xinit);
}
void mass_resetF(t_mass *x)
{
x->force=0;
}
void mass_dX(t_mass *x, t_float posX)
{
x->dX += posX;
}
void mass_setX(t_mass *x, t_float posX)
{
x->pos_old_2 = posX; // clear history for stability (instability) problem
x->pos_old_1 = posX;
x->force=0;
outlet_float(x->position_new, posX);
}
void mass_loadbang(t_mass *x)
{
outlet_float(x->position_new, x->Xinit);
}
void mass_set_mass(t_mass *x, t_float mass)
{
x->mass=mass;
}
static void mass_free(t_mass *x)
{
pd_unbind(&x->x_obj.ob_pd, x->x_sym);
}
void *mass_new(t_symbol *s, t_floatarg M, t_floatarg X)
{
t_mass *x = (t_mass *)pd_new(mass_class);
x->x_sym = s;
pd_bind(&x->x_obj.ob_pd, s);
x->position_new=outlet_new(&x->x_obj, 0);
x->force_out=outlet_new(&x->x_obj, 0);
x->vitesse_out=outlet_new(&x->x_obj, 0);
x->Xinit=X;
x->pos_old_1 = X;
x->pos_old_2 = X;
x->force=0;
x->mass=M;
x->minX = -100000;
x->maxX = 100000;
if (x->mass<=0) x->mass=1;
makeseed();
return (void *)x;
}
void mass_setup(void)
{
mass_class = class_new(gensym("mass"),
(t_newmethod)mass_new,
(t_method)mass_free,
sizeof(t_mass),
CLASS_DEFAULT, A_DEFSYM, A_DEFFLOAT, A_DEFFLOAT,0);
class_addcreator((t_newmethod)mass_new, gensym("masse"), A_DEFSYM, A_DEFFLOAT, A_DEFFLOAT,0);
class_addfloat(mass_class, mass_float);
class_addbang(mass_class, mass_bang);
class_addmethod(mass_class, (t_method)mass_set_mass, gensym("setM"), A_DEFFLOAT, 0);
class_addmethod(mass_class, (t_method)mass_setX, gensym("setX"), A_DEFFLOAT, 0);
class_addmethod(mass_class, (t_method)mass_dX, gensym("dX"), A_DEFFLOAT, 0);
class_addmethod(mass_class, (t_method)mass_reset, gensym("reset"), 0);
class_addmethod(mass_class, (t_method)mass_resetF, gensym("resetF"), 0);
class_addmethod(mass_class, (t_method)mass_minX, gensym("setXmin"), A_DEFFLOAT, 0);
class_addmethod(mass_class, (t_method)mass_maxX, gensym("setXmax"), A_DEFFLOAT, 0);
class_addmethod(mass_class, (t_method)mass_loadbang, gensym("loadbang"), 0);
}
[/code]
[code]
#include "m_pd.h"
#include "math.h"
static t_class *linkKD_class;
typedef struct _linkKD {
t_object x_obj;
t_float raideur, viscosite, D2, longueur, distance_old, position1, position2, position_old1, position_old2;
t_outlet *force1;
t_outlet *force2;
t_float Lmin, Lmax;
t_symbol *x_sym; // receive
} t_linkKD;
void linkKD_float(t_linkKD *x, t_floatarg f1)
{
x->position1 = f1;
}
void linkKD_bang(t_linkKD *x)
{
t_float force1, force2, distance;
distance = (x->position2 - x->position1);
//distance = abs(x->position2 - x->position1);
if (distance<0) distance = -distance;
force1 = x->raideur*(distance-(x->longueur)) + x->viscosite*(distance - x->distance_old) ;
x->distance_old = distance;
if (distance > x->Lmax) force1=0;
if (distance < x->Lmin) force1=0;
if (distance != 0)
{
force1 = force1 * (x->position2 - x->position1) / distance;
}
force2 = -force1 + (x->position_old2 - x->position2)*x->D2;
force1 += (x->position_old1 - x->position1)*x->D2;
// masse damping
outlet_float(x->force1, force1);
outlet_float(x->force2, force2);
x->position_old1 = x->position1;
x->position_old2 = x->position2;
}
void linkKD_reset(t_linkKD *x)
{
x->position1 = 0;
x->position2 = 0;
x->position_old1 = 0;
x->position_old2 = 0;
x->distance_old = x->longueur;
}
void linkKD_resetF(t_linkKD *x)
{
x->position_old1 = x->position1;
x->position_old2 = x->position2;
x->distance_old = x->longueur;
}
void linkKD_resetl(t_linkKD *x)
{
x->longueur = (x->position1 - x->position2);
}
void linkKD_setL(t_linkKD *x, t_float L)
{
x->longueur = L;
}
void linkKD_setK(t_linkKD *x, t_float K)
{
x->raideur = K;
}
void linkKD_setD(t_linkKD *x, t_float D)
{
x->viscosite = D;
}
void linkKD_setD2(t_linkKD *x, t_float D2)
{
x->D2 = D2;
}
void linkKD_Lmin(t_linkKD *x, t_float Lmin)
{
x->Lmin = Lmin;
}
void linkKD_Lmax(t_linkKD *x, t_float Lmax)
{
x->Lmax = Lmax;
}
static void linkKD_free(t_linkKD *x)
{
pd_unbind(&x->x_obj.ob_pd, x->x_sym);
}
void *linkKD_new(t_symbol *s, t_floatarg L, t_floatarg K, t_floatarg D, t_floatarg D2 )
{
t_linkKD *x = (t_linkKD *)pd_new(linkKD_class);
x->x_sym = s;
pd_bind(&x->x_obj.ob_pd, s);
floatinlet_new(&x->x_obj, &x->position2);
x->force1=outlet_new(&x->x_obj, 0);
x->force2=outlet_new(&x->x_obj, 0);
x->position1 = 0;
x->position2 = 0;
x->raideur=K;
x->viscosite=D;
x->D2=D2;
x->Lmin= 0;
x->Lmax= 10000;
x->longueur=L;
return (void *)x;
}
void link_setup(void)
{
linkKD_class = class_new(gensym("link"),
(t_newmethod)linkKD_new,
(t_method)linkKD_free,
sizeof(t_linkKD),
CLASS_DEFAULT, A_DEFSYM, A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0);
class_addcreator((t_newmethod)linkKD_new, gensym("lia"), A_DEFSYM, A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, A_DEFFLOAT, 0);
class_addfloat(linkKD_class, linkKD_float);
class_addbang(linkKD_class, linkKD_bang);
class_addmethod(linkKD_class, (t_method)linkKD_reset, gensym("reset"), 0);
class_addmethod(linkKD_class, (t_method)linkKD_resetl, gensym("resetL"), 0);
class_addmethod(linkKD_class, (t_method)linkKD_resetF, gensym("resetF"), 0);
class_addmethod(linkKD_class, (t_method)linkKD_setD, gensym("setD"), A_DEFFLOAT, 0);
class_addmethod(linkKD_class, (t_method)linkKD_setD2, gensym("setD2"), A_DEFFLOAT, 0);
class_addmethod(linkKD_class, (t_method)linkKD_setK, gensym("setK"), A_DEFFLOAT, 0);
class_addmethod(linkKD_class, (t_method)linkKD_setL, gensym("setL"), A_DEFFLOAT, 0);
class_addmethod(linkKD_class, (t_method)linkKD_Lmin, gensym("setLmin"), A_DEFFLOAT, 0);
class_addmethod(linkKD_class, (t_method)linkKD_Lmax, gensym("setLmax"), A_DEFFLOAT, 0);
}
[/code]
i'm not running in to those instability issues, but it's difficult to say why from your video. you may want to set the "resetquat" attribute to align with the rotation that is required by the barslide constraint (this will eliminate the initial rotation as the bodies are forced to align themselves to the constraint). if you are still getting those rotation errors it should be easy enough to zero them out every frame.