|
|
|
|
@@ -87,6 +87,7 @@ variables on the UI for now
|
|
|
|
|
#include "BIF_editdeform.h"
|
|
|
|
|
#include "BIF_graphics.h"
|
|
|
|
|
#include "PIL_time.h"
|
|
|
|
|
#include "ONL_opennl.h"
|
|
|
|
|
|
|
|
|
|
/* callbacks for errors and interrupts and some goo */
|
|
|
|
|
static int (*SB_localInterruptCallBack)(void) = NULL;
|
|
|
|
|
@@ -97,7 +98,7 @@ static int (*SB_localInterruptCallBack)(void) = NULL;
|
|
|
|
|
|
|
|
|
|
typedef struct BodySpring {
|
|
|
|
|
int v1, v2;
|
|
|
|
|
float len, strength, cf;
|
|
|
|
|
float len, strength, cf,load;
|
|
|
|
|
float ext_force[3]; /* edges colliding and sailing */
|
|
|
|
|
short order;
|
|
|
|
|
short flag;
|
|
|
|
|
@@ -120,6 +121,11 @@ typedef struct SBScratch {
|
|
|
|
|
float aabbmin[3],aabbmax[3];
|
|
|
|
|
}SBScratch;
|
|
|
|
|
|
|
|
|
|
#define NLF_BUILD 1
|
|
|
|
|
#define NLF_SOLVE 2
|
|
|
|
|
|
|
|
|
|
#define MID_PRESERVE 1
|
|
|
|
|
|
|
|
|
|
#define SOFTGOALSNAP 0.999f
|
|
|
|
|
/* if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
|
|
|
|
|
removes *unnecessary* stiffnes from ODE system
|
|
|
|
|
@@ -585,6 +591,7 @@ static void add_mesh_quad_diag_springs(Object *ob)
|
|
|
|
|
|
|
|
|
|
if (ob->soft){
|
|
|
|
|
int nofquads;
|
|
|
|
|
float s_shear = ob->soft->shearstiff*ob->soft->shearstiff;
|
|
|
|
|
|
|
|
|
|
nofquads = count_mesh_quads(me);
|
|
|
|
|
if (nofquads) {
|
|
|
|
|
@@ -605,12 +612,12 @@ static void add_mesh_quad_diag_springs(Object *ob)
|
|
|
|
|
if(mface->v4) {
|
|
|
|
|
bs->v1= mface->v1;
|
|
|
|
|
bs->v2= mface->v3;
|
|
|
|
|
bs->strength= 1.0;
|
|
|
|
|
bs->strength= s_shear;
|
|
|
|
|
bs->order =2;
|
|
|
|
|
bs++;
|
|
|
|
|
bs->v1= mface->v2;
|
|
|
|
|
bs->v2= mface->v4;
|
|
|
|
|
bs->strength= 1.0;
|
|
|
|
|
bs->strength= s_shear;
|
|
|
|
|
bs->order =2;
|
|
|
|
|
bs++;
|
|
|
|
|
|
|
|
|
|
@@ -692,6 +699,7 @@ static void add_2nd_order_springs(Object *ob,float stiffness)
|
|
|
|
|
{
|
|
|
|
|
int counter = 0;
|
|
|
|
|
BodySpring *bs_new;
|
|
|
|
|
stiffness *=stiffness;
|
|
|
|
|
|
|
|
|
|
add_2nd_order_roller(ob,stiffness,&counter,0); /* counting */
|
|
|
|
|
if (counter) {
|
|
|
|
|
@@ -1450,6 +1458,7 @@ int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
void scan_for_ext_spring_forces(Object *ob,float timenow)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb = ob->soft;
|
|
|
|
|
@@ -2095,8 +2104,110 @@ static int sb_deflect_face(Object *ob,float *actpos,float *facenormal,float *for
|
|
|
|
|
return(deflected);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
static void dfdx_spring(int ia, int ic, int op, float dir[3],float L,float len,float factor)
|
|
|
|
|
{
|
|
|
|
|
float m,delta_ij;
|
|
|
|
|
int i ,j;
|
|
|
|
|
if (L < len){
|
|
|
|
|
for(i=0;i<3;i++)
|
|
|
|
|
for(j=0;j<3;j++){
|
|
|
|
|
delta_ij = (i==j ? (1.0f): (0.0f));
|
|
|
|
|
m=factor*(dir[i]*dir[j] + (1-L/len)*(delta_ij - dir[i]*dir[j]));
|
|
|
|
|
nlMatrixAdd(ia+i,op+ic+j,m);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
else{
|
|
|
|
|
for(i=0;i<3;i++)
|
|
|
|
|
for(j=0;j<3;j++){
|
|
|
|
|
m=factor*dir[i]*dir[j];
|
|
|
|
|
nlMatrixAdd(ia+i,op+ic+j,m);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
|
|
|
|
|
static void dfdx_goal(int ia, int ic, int op, float factor)
|
|
|
|
|
{
|
|
|
|
|
int i;
|
|
|
|
|
for(i=0;i<3;i++) nlMatrixAdd(ia+i,op+ic+i,factor);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
static void dfdv_goal(int ia, int ic,float factor)
|
|
|
|
|
{
|
|
|
|
|
int i;
|
|
|
|
|
for(i=0;i<3;i++) nlMatrixAdd(ia+i,ic+i,factor);
|
|
|
|
|
}
|
|
|
|
|
static void sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there */
|
|
|
|
|
BodyPoint *bp1,*bp2;
|
|
|
|
|
|
|
|
|
|
float dir[3],dvel[3];
|
|
|
|
|
float distance,forcefactor,kd,absvel,projvel;
|
|
|
|
|
int ia,ic;
|
|
|
|
|
/* prepare depending on which side of the spring we are on */
|
|
|
|
|
if (bpi == bs->v1){
|
|
|
|
|
bp1 = &sb->bpoint[bs->v1];
|
|
|
|
|
bp2 = &sb->bpoint[bs->v2];
|
|
|
|
|
ia =3*bs->v1;
|
|
|
|
|
ic =3*bs->v2;
|
|
|
|
|
}
|
|
|
|
|
else if (bpi == bs->v2){
|
|
|
|
|
bp1 = &sb->bpoint[bs->v2];
|
|
|
|
|
bp2 = &sb->bpoint[bs->v1];
|
|
|
|
|
ia =3*bs->v2;
|
|
|
|
|
ic =3*bs->v1;
|
|
|
|
|
}
|
|
|
|
|
else{
|
|
|
|
|
/* TODO make this debug option */
|
|
|
|
|
/**/
|
|
|
|
|
printf("bodypoint <bpi> is not attached to spring <*bs> --> sb_spring_force()\n");
|
|
|
|
|
return;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* do bp1 <--> bp2 elastic */
|
|
|
|
|
VecSubf(dir,bp1->pos,bp2->pos);
|
|
|
|
|
distance = Normalize(dir);
|
|
|
|
|
if (bs->len < distance)
|
|
|
|
|
iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
|
|
|
|
|
else
|
|
|
|
|
iks = 1.0f/(1.0f-sb->inpush)-1.0f ;/* inner spring constants function */
|
|
|
|
|
|
|
|
|
|
if(bs->len > 0.0f) /* check for degenerated springs */
|
|
|
|
|
forcefactor = iks/bs->len;
|
|
|
|
|
else
|
|
|
|
|
forcefactor = iks;
|
|
|
|
|
forcefactor *= bs->strength;
|
|
|
|
|
Vec3PlusStVec(bp1->force,(bs->len - distance)*forcefactor,dir);
|
|
|
|
|
|
|
|
|
|
/* do bp1 <--> bp2 viscous */
|
|
|
|
|
VecSubf(dvel,bp1->vec,bp2->vec);
|
|
|
|
|
kd = sb->infrict * sb_fric_force_scale(ob);
|
|
|
|
|
absvel = Normalize(dvel);
|
|
|
|
|
projvel = Inpf(dir,dvel);
|
|
|
|
|
kd *= absvel * projvel;
|
|
|
|
|
Vec3PlusStVec(bp1->force,-kd,dir);
|
|
|
|
|
|
|
|
|
|
/* do jacobian stuff if needed */
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int op =3*sb->totpoint;
|
|
|
|
|
float mvel = -forcetime*kd;
|
|
|
|
|
float mpos = -forcetime*forcefactor;
|
|
|
|
|
/* depending on my pos */
|
|
|
|
|
dfdx_spring(ia,ia,op,dir,bs->len,distance,-mpos);
|
|
|
|
|
/* depending on my vel */
|
|
|
|
|
dfdv_goal(ia,ia,mvel); // well that ignores geometie
|
|
|
|
|
if(bp2->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */
|
|
|
|
|
/* depending on other pos */
|
|
|
|
|
dfdx_spring(ia,ic,op,dir,bs->len,distance,mpos);
|
|
|
|
|
/* depending on other vel */
|
|
|
|
|
dfdv_goal(ia,ia,-mvel); // well that ignores geometie
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags)
|
|
|
|
|
{
|
|
|
|
|
/* rule we never alter free variables :bp->vec bp->pos in here !
|
|
|
|
|
* this will ruin adaptive stepsize AKA heun! (BM)
|
|
|
|
|
@@ -2106,10 +2217,20 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
BodyPoint *bproot;
|
|
|
|
|
BodySpring *bs;
|
|
|
|
|
ListBase *do_effector;
|
|
|
|
|
float iks, ks, kd, gravity, actspringlen, forcefactor, sd[3];
|
|
|
|
|
float iks, ks, kd, gravity;
|
|
|
|
|
float fieldfactor = 1000.0f, windfactor = 250.0f;
|
|
|
|
|
float tune = sb->ballstiff;
|
|
|
|
|
int a, b, do_deflector,do_selfcollision,do_springcollision,do_aero;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
NLboolean success;
|
|
|
|
|
|
|
|
|
|
if(nl_flags){
|
|
|
|
|
nlBegin(NL_SYSTEM);
|
|
|
|
|
nlBegin(NL_MATRIX);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
gravity = sb->grav * sb_grav_force_scale(ob);
|
|
|
|
|
@@ -2131,7 +2252,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
float defforce[3];
|
|
|
|
|
do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/*
|
|
|
|
|
if (do_selfcollision ){
|
|
|
|
|
float ce[3];
|
|
|
|
|
VecMidf(ce,sb->scratch->aabbmax,sb->scratch->aabbmin);
|
|
|
|
|
@@ -2139,10 +2260,33 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
bp->octantflag = set_octant_flags(ce,bp->pos,bp->colball);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
*/
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
/* clear forces accumulator */
|
|
|
|
|
bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
int op =3*sb->totpoint;
|
|
|
|
|
/* dF/dV = v */
|
|
|
|
|
|
|
|
|
|
nlMatrixAdd(op+ia,ia,-forcetime);
|
|
|
|
|
nlMatrixAdd(op+ia+1,ia+1,-forcetime);
|
|
|
|
|
nlMatrixAdd(op+ia+2,ia+2,-forcetime);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* we need [unit - h* dF/dy]^-1 */
|
|
|
|
|
|
|
|
|
|
nlMatrixAdd(ia,ia,1);
|
|
|
|
|
nlMatrixAdd(ia+1,ia+1,1);
|
|
|
|
|
nlMatrixAdd(ia+2,ia+2,1);
|
|
|
|
|
|
|
|
|
|
nlMatrixAdd(op+ia,op+ia,1);
|
|
|
|
|
nlMatrixAdd(op+ia+1,op+ia+1,1);
|
|
|
|
|
nlMatrixAdd(op+ia+2,op+ia+2,1);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* naive ball self collision */
|
|
|
|
|
/* needs to be done if goal snaps or not */
|
|
|
|
|
@@ -2156,7 +2300,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
|
|
|
|
|
for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) {
|
|
|
|
|
|
|
|
|
|
if ((bp->octantflag & obp->octantflag) == 0) continue;
|
|
|
|
|
//if ((bp->octantflag & obp->octantflag) == 0) continue;
|
|
|
|
|
|
|
|
|
|
compare = (obp->colball + bp->colball);
|
|
|
|
|
VecSubf(def, bp->pos, obp->pos);
|
|
|
|
|
@@ -2182,8 +2326,33 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
VecSubf(dvel,velcenter,bp->vec);
|
|
|
|
|
VecMulf(dvel,sb->nodemass);
|
|
|
|
|
|
|
|
|
|
Vec3PlusStVec(bp->force,sb->balldamp,dvel);
|
|
|
|
|
Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
|
|
|
|
|
Vec3PlusStVec(bp->force,sb->balldamp,dvel);
|
|
|
|
|
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
int ic =3*(sb->totpoint-c);
|
|
|
|
|
int op =3*sb->totpoint;
|
|
|
|
|
float mvel = forcetime*sb->nodemass*sb->balldamp;
|
|
|
|
|
float mpos = forcetime*tune*(1.0f-sb->balldamp);
|
|
|
|
|
/*some quick and dirty entries to the jacobian*/
|
|
|
|
|
dfdx_goal(ia,ia,op,mpos);
|
|
|
|
|
dfdv_goal(ia,ia,mvel);
|
|
|
|
|
/* exploit force(a,b) == -force(b,a) part1/2 */
|
|
|
|
|
dfdx_goal(ic,ic,op,mpos);
|
|
|
|
|
dfdv_goal(ic,ic,mvel);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*TODO sit down an X-out the true jacobian entries*/
|
|
|
|
|
/*well does not make to much sense because the eigenvalues
|
|
|
|
|
of the jacobian go negative; and negative eigenvalues
|
|
|
|
|
on a complex iterative system z(n+1)=A * z(n)
|
|
|
|
|
give imaginary roots in the charcateristic polynom
|
|
|
|
|
--> solutions that to z(t)=u(t)* exp ( i omega t) --> oscilations we don't want here
|
|
|
|
|
where u(t) is a unknown amplitude function (worst case rising fast)
|
|
|
|
|
*/
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* exploit force(a,b) == -force(b,a) part2/2 */
|
|
|
|
|
VecSubf(dvel,velcenter,obp->vec);
|
|
|
|
|
VecMulf(dvel,sb->nodemass);
|
|
|
|
|
@@ -2191,6 +2360,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
Vec3PlusStVec(obp->force,sb->balldamp,dvel);
|
|
|
|
|
Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
@@ -2201,23 +2371,39 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
float auxvect[3];
|
|
|
|
|
float velgoal[3];
|
|
|
|
|
float absvel =0, projvel= 0;
|
|
|
|
|
|
|
|
|
|
/* do goal stuff */
|
|
|
|
|
if(ob->softflag & OB_SB_GOAL) {
|
|
|
|
|
/* true elastic goal */
|
|
|
|
|
VecSubf(auxvect,bp->origT,bp->pos);
|
|
|
|
|
VecSubf(auxvect,bp->pos,bp->origT);
|
|
|
|
|
ks = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
|
|
|
|
|
bp->force[0]+= ks*(auxvect[0]);
|
|
|
|
|
bp->force[1]+= ks*(auxvect[1]);
|
|
|
|
|
bp->force[2]+= ks*(auxvect[2]);
|
|
|
|
|
bp->force[0]+= -ks*(auxvect[0]);
|
|
|
|
|
bp->force[1]+= -ks*(auxvect[1]);
|
|
|
|
|
bp->force[2]+= -ks*(auxvect[2]);
|
|
|
|
|
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
int op =3*(sb->totpoint);
|
|
|
|
|
/* depending on my pos */
|
|
|
|
|
dfdx_goal(ia,ia,op,ks*forcetime);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* calulate damping forces generated by goals*/
|
|
|
|
|
VecSubf(velgoal,bp->origS, bp->origE);
|
|
|
|
|
kd = sb->goalfrict * sb_fric_force_scale(ob) ;
|
|
|
|
|
VecAddf(auxvect,velgoal,bp->vec);
|
|
|
|
|
|
|
|
|
|
if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */
|
|
|
|
|
bp->force[0]-= kd * (velgoal[0] + bp->vec[0]);
|
|
|
|
|
bp->force[1]-= kd * (velgoal[1] + bp->vec[1]);
|
|
|
|
|
bp->force[2]-= kd * (velgoal[2] + bp->vec[2]);
|
|
|
|
|
bp->force[0]-= kd * (auxvect[0]);
|
|
|
|
|
bp->force[1]-= kd * (auxvect[1]);
|
|
|
|
|
bp->force[2]-= kd * (auxvect[2]);
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
Normalize(auxvect);
|
|
|
|
|
/* depending on my vel */
|
|
|
|
|
dfdv_goal(ia,ia,kd*forcetime);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
|
|
|
|
|
@@ -2230,6 +2416,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
|
|
|
|
|
/* gravitation */
|
|
|
|
|
bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
|
|
|
|
|
//bp->force[1]-= gravity*sb->nodemass; /* individual mass of node here */
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* particle field & vortex */
|
|
|
|
|
@@ -2260,6 +2447,16 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
bp->force[1]-= bp->vec[1]*kd;
|
|
|
|
|
bp->force[2]-= bp->vec[2]*kd;
|
|
|
|
|
/* friction in media done */
|
|
|
|
|
if(nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
int op =3*sb->totpoint;
|
|
|
|
|
/* da/dv = */
|
|
|
|
|
|
|
|
|
|
nlMatrixAdd(ia,ia,forcetime*kd);
|
|
|
|
|
nlMatrixAdd(ia+1,ia+1,forcetime*kd);
|
|
|
|
|
nlMatrixAdd(ia+2,ia+2,forcetime*kd);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
/* +++cached collision targets */
|
|
|
|
|
bp->choke = 0.0f;
|
|
|
|
|
@@ -2269,7 +2466,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
kd = 1.0f;
|
|
|
|
|
|
|
|
|
|
if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
|
|
|
|
|
if (intrusion < 0.0f){
|
|
|
|
|
if ((!nl_flags)&&(intrusion < 0.0f)){
|
|
|
|
|
/*bjornmose: uugh.. what an evil hack
|
|
|
|
|
violation of the 'don't touch bp->pos in here' rule
|
|
|
|
|
but works nice, like this-->
|
|
|
|
|
@@ -2287,8 +2484,15 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
VECSUB(cfforce,bp->vec,vel);
|
|
|
|
|
Vec3PlusStVec(bp->force,-cf*50.0f,cfforce);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
Vec3PlusStVec(bp->force,kd,defforce);
|
|
|
|
|
if (nl_flags & NLF_BUILD){
|
|
|
|
|
int ia =3*(sb->totpoint-a);
|
|
|
|
|
int op =3*sb->totpoint;
|
|
|
|
|
//dfdx_goal(ia,ia,op,mpos); // don't do unless you know
|
|
|
|
|
//dfdv_goal(ia,ia,-cf);
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
@@ -2305,48 +2509,8 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
bp->choke = bs->cf;
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (( (sb->totpoint-a) == bs->v1) ){
|
|
|
|
|
VecSubf(sd,(bproot+bs->v2)->pos, bp->pos);
|
|
|
|
|
actspringlen=Normalize(sd);
|
|
|
|
|
|
|
|
|
|
/* friction stuff V1 */
|
|
|
|
|
VecSubf(velgoal,bp->vec,(bproot+bs->v2)->vec);
|
|
|
|
|
kd = sb->infrict * sb_fric_force_scale(ob);
|
|
|
|
|
absvel = Normalize(velgoal);
|
|
|
|
|
projvel = ABS(Inpf(sd,velgoal));
|
|
|
|
|
kd *= absvel * projvel;
|
|
|
|
|
Vec3PlusStVec(bp->force,-kd,velgoal);
|
|
|
|
|
|
|
|
|
|
if(bs->len > 0.0f) /* check for degenerated springs */
|
|
|
|
|
forcefactor = (bs->len - actspringlen)/bs->len * iks;
|
|
|
|
|
else
|
|
|
|
|
forcefactor = actspringlen * iks;
|
|
|
|
|
forcefactor *= bs->strength;
|
|
|
|
|
|
|
|
|
|
Vec3PlusStVec(bp->force,-forcefactor,sd);
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (( (sb->totpoint-a) == bs->v2) ){
|
|
|
|
|
VecSubf(sd,bp->pos,(bproot+bs->v1)->pos);
|
|
|
|
|
actspringlen=Normalize(sd);
|
|
|
|
|
|
|
|
|
|
/* friction stuff V2 */
|
|
|
|
|
VecSubf(velgoal,bp->vec,(bproot+bs->v1)->vec);
|
|
|
|
|
kd = sb->infrict * sb_fric_force_scale(ob);
|
|
|
|
|
absvel = Normalize(velgoal);
|
|
|
|
|
projvel = ABS(Inpf(sd,velgoal));
|
|
|
|
|
kd *= absvel * projvel;
|
|
|
|
|
Vec3PlusStVec(bp->force,-kd,velgoal);
|
|
|
|
|
|
|
|
|
|
if(bs->len > 0.0)
|
|
|
|
|
forcefactor = (bs->len - actspringlen)/bs->len * iks;
|
|
|
|
|
else
|
|
|
|
|
forcefactor = actspringlen * iks;
|
|
|
|
|
forcefactor *= bs->strength;
|
|
|
|
|
Vec3PlusStVec(bp->force,+forcefactor,sd);
|
|
|
|
|
}
|
|
|
|
|
// sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
|
|
|
|
|
sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,nl_flags);
|
|
|
|
|
}/* loop springs */
|
|
|
|
|
}/* existing spring list */
|
|
|
|
|
}/*any edges*/
|
|
|
|
|
@@ -2356,14 +2520,206 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* finally add forces caused by face collision */
|
|
|
|
|
if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
|
|
|
|
|
if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
|
|
|
|
|
|
|
|
|
|
/* finish matrix and solve */
|
|
|
|
|
if(nl_flags & NLF_SOLVE){
|
|
|
|
|
//double sct,sst=PIL_check_seconds_timer();
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
int iv =3*(sb->totpoint-a);
|
|
|
|
|
int ip =3*(2*sb->totpoint-a);
|
|
|
|
|
int n;
|
|
|
|
|
for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);}
|
|
|
|
|
for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);}
|
|
|
|
|
}
|
|
|
|
|
nlEnd(NL_MATRIX);
|
|
|
|
|
nlEnd(NL_SYSTEM);
|
|
|
|
|
|
|
|
|
|
if ((G.rt >0) && (nl_flags & NLF_BUILD))
|
|
|
|
|
{
|
|
|
|
|
printf("####MEE#####\n");
|
|
|
|
|
nlPrintMatrix();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
success= nlSolveAdvanced(NULL, 1);
|
|
|
|
|
|
|
|
|
|
// nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */
|
|
|
|
|
if(success){
|
|
|
|
|
float f;
|
|
|
|
|
int index =0;
|
|
|
|
|
/* for debug purpose .. anyhow cropping B vector looks like working */
|
|
|
|
|
if (G.rt >0)
|
|
|
|
|
for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
printf("(%f ",f);index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
printf("%f ",f);index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
printf("%f)",f);index++;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
index =0;
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdv[0] = f; index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdv[1] = f; index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdv[2] = f; index++;
|
|
|
|
|
}
|
|
|
|
|
/*
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdx[0] = f; index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdx[1] = f; index++;
|
|
|
|
|
f=nlGetVariable(0,index);
|
|
|
|
|
bp->impdx[2] = f; index++;
|
|
|
|
|
}
|
|
|
|
|
*/
|
|
|
|
|
}
|
|
|
|
|
else{
|
|
|
|
|
printf("Matrix inversion failed \n");
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
VECCOPY(bp->impdv,bp->force);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//sct=PIL_check_seconds_timer();
|
|
|
|
|
//if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name);
|
|
|
|
|
}
|
|
|
|
|
/* cleanup */
|
|
|
|
|
|
|
|
|
|
if(do_effector) pdEndEffectors(do_effector);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
static void softbody_apply_fake_implicit_forces(Object *ob, float forcetime, float *err, float *err_ref_pos,float *err_ref_vel)
|
|
|
|
|
{
|
|
|
|
|
/* time evolution */
|
|
|
|
|
/* do semi implicit euler */
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there */
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
float *perp,*perv;
|
|
|
|
|
float erp[3],erv[3],dx[3],dv[3],aabbmin[3],aabbmax[3],cm[3]={0.0f,0.0f,0.0f};
|
|
|
|
|
float timeovermass;
|
|
|
|
|
float maxerrpos= 0.0f,maxerrvel = 0.0f,maxerrvel2 = 0.0f;
|
|
|
|
|
int a,fuzzy=0;
|
|
|
|
|
|
|
|
|
|
static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err)
|
|
|
|
|
forcetime *= sb_time_scale(ob);
|
|
|
|
|
|
|
|
|
|
aabbmin[0]=aabbmin[1]=aabbmin[2] = 1e20f;
|
|
|
|
|
aabbmax[0]=aabbmax[1]=aabbmax[2] = -1e20f;
|
|
|
|
|
|
|
|
|
|
/* claim a minimum mass for vertex */
|
|
|
|
|
if (sb->nodemass > 0.009999f) timeovermass = forcetime/sb->nodemass;
|
|
|
|
|
else timeovermass = forcetime/0.009999f;
|
|
|
|
|
/* step along error reference vector */
|
|
|
|
|
perp =err_ref_pos;
|
|
|
|
|
perv =err_ref_vel;
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
/* step along error reference vector */
|
|
|
|
|
if(perp) {VECCOPY(erp,perp);perp +=3;}
|
|
|
|
|
if(perv) {VECCOPY(erv,perv);perv +=3;}
|
|
|
|
|
if(bp->goal < SOFTGOALSNAP){
|
|
|
|
|
|
|
|
|
|
/* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces + more forces*/
|
|
|
|
|
/* the ( ... )' operator denotes derivate respective time */
|
|
|
|
|
/* the "implicit" euler step for velocity then becomes */
|
|
|
|
|
/* v(t + dt) = v(t) + O * X(t) * dt */
|
|
|
|
|
/* O = [1-dt*dA/d(X)]^1 */
|
|
|
|
|
/* with X = (a[n],v[n]) */
|
|
|
|
|
/* da/dv | da/dx */
|
|
|
|
|
/* dA/d(X) = ( ------------- )*/
|
|
|
|
|
/* dv/dv | dv/dx */
|
|
|
|
|
/* note!! we did not solve any implicit system but looking forward more or less smart
|
|
|
|
|
what a implicit solution may look like by means of taylor expansion */
|
|
|
|
|
VECCOPY(dx,bp->vec);
|
|
|
|
|
|
|
|
|
|
bp->force[0]= timeovermass * bp->impdv[0]; /* individual mass of node here */
|
|
|
|
|
bp->force[1]= timeovermass * bp->impdv[1];
|
|
|
|
|
bp->force[2]= timeovermass * bp->impdv[2];
|
|
|
|
|
VECCOPY(dv,bp->force);
|
|
|
|
|
if(perv){
|
|
|
|
|
maxerrvel = MAX2(maxerrvel,ABS(dv[0] - erv[0]));
|
|
|
|
|
maxerrvel = MAX2(maxerrvel,ABS(dv[1] - erv[1]));
|
|
|
|
|
maxerrvel = MAX2(maxerrvel,ABS(dv[2] - erv[2]));
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
maxerrvel2 = MAX2(maxerrvel2,ABS(dv[0]));
|
|
|
|
|
maxerrvel2 = MAX2(maxerrvel2,ABS(dv[1]));
|
|
|
|
|
maxerrvel2 = MAX2(maxerrvel2,ABS(dv[2]));
|
|
|
|
|
|
|
|
|
|
VECADD(bp->vec, bp->vec, dv);
|
|
|
|
|
|
|
|
|
|
/* so here is (x)'= v(elocity) */
|
|
|
|
|
/* the euler step for location then becomes */
|
|
|
|
|
/* x(t + dt) = x(t) + v(t) * dt */
|
|
|
|
|
|
|
|
|
|
// VECCOPY(dx,bp->vec);
|
|
|
|
|
|
|
|
|
|
dx[0]*=forcetime;
|
|
|
|
|
dx[1]*=forcetime;
|
|
|
|
|
dx[2]*=forcetime;
|
|
|
|
|
|
|
|
|
|
/* bp->choke is set when we need to pull a vertex or edge out of the collider.
|
|
|
|
|
the collider object signals to get out by pushing hard. on the other hand
|
|
|
|
|
we don't want to end up in deep space so we add some <viscosity>
|
|
|
|
|
to balance that out */
|
|
|
|
|
if (bp->choke > 0.0f){
|
|
|
|
|
dx[0]*=(1.0f - bp->choke);
|
|
|
|
|
dx[1]*=(1.0f - bp->choke);
|
|
|
|
|
dx[2]*=(1.0f - bp->choke);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* should be possible to get more out of the matrix inversion
|
|
|
|
|
but not verified yet
|
|
|
|
|
dx[0]*=forcetime*forcetime*bp->impdx[0];
|
|
|
|
|
dx[1]*=forcetime*forcetime*bp->impdx[1];
|
|
|
|
|
dx[2]*=forcetime*forcetime*bp->impdx[2];
|
|
|
|
|
*/
|
|
|
|
|
VECADD(bp->pos, bp->pos, dx);
|
|
|
|
|
if (perp){
|
|
|
|
|
maxerrpos = MAX2(maxerrpos,ABS(bp->pos[0]-erp[0]));
|
|
|
|
|
maxerrpos = MAX2(maxerrpos,ABS(bp->pos[1]-erp[1]));
|
|
|
|
|
maxerrpos = MAX2(maxerrpos,ABS(bp->pos[2]-erp[2]));
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}/*snap*/
|
|
|
|
|
/* so while we are looping BPs anyway do statistics on the fly */
|
|
|
|
|
aabbmin[0] = MIN2(aabbmin[0],bp->pos[0]);
|
|
|
|
|
aabbmin[1] = MIN2(aabbmin[1],bp->pos[1]);
|
|
|
|
|
aabbmin[2] = MIN2(aabbmin[2],bp->pos[2]);
|
|
|
|
|
aabbmax[0] = MAX2(aabbmax[0],bp->pos[0]);
|
|
|
|
|
aabbmax[1] = MAX2(aabbmax[1],bp->pos[1]);
|
|
|
|
|
aabbmax[2] = MAX2(aabbmax[2],bp->pos[2]);
|
|
|
|
|
if (bp->flag & SBF_DOFUZZY) fuzzy =1;
|
|
|
|
|
} /*for*/
|
|
|
|
|
|
|
|
|
|
if (sb->totpoint) VecMulf(cm,1.0f/sb->totpoint);
|
|
|
|
|
if (sb->scratch){
|
|
|
|
|
VECCOPY(sb->scratch->aabbmin,aabbmin);
|
|
|
|
|
VECCOPY(sb->scratch->aabbmax,aabbmax);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (err){ /* so step size will be controlled by biggest difference in slope */
|
|
|
|
|
if (sb->solverflags & SBSO_OLDERR)
|
|
|
|
|
*err = MAX2(maxerrpos,maxerrvel);
|
|
|
|
|
else
|
|
|
|
|
if (perp) *err = maxerrpos;
|
|
|
|
|
else *err = maxerrvel2;
|
|
|
|
|
//printf("EP %f EV %f \n",maxerrpos,maxerrvel);
|
|
|
|
|
if (fuzzy){
|
|
|
|
|
*err /= sb->fuzzyness;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags)
|
|
|
|
|
{
|
|
|
|
|
/* time evolution */
|
|
|
|
|
/* actually does an explicit euler step mode == 0 */
|
|
|
|
|
@@ -2386,6 +2742,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
if(bp->goal < SOFTGOALSNAP){
|
|
|
|
|
if(mid_flags & MID_PRESERVE) VECCOPY(dx,bp->vec);
|
|
|
|
|
|
|
|
|
|
/* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces + more forces*/
|
|
|
|
|
/* the ( ... )' operator denotes derivate respective time */
|
|
|
|
|
@@ -2417,8 +2774,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *
|
|
|
|
|
/* so here is (x)'= v(elocity) */
|
|
|
|
|
/* the euler step for location then becomes */
|
|
|
|
|
/* x(t + dt) = x(t) + v(t) * dt */
|
|
|
|
|
|
|
|
|
|
VECCOPY(dx,bp->vec);
|
|
|
|
|
if(!(mid_flags & MID_PRESERVE)) VECCOPY(dx,bp->vec);
|
|
|
|
|
dx[0]*=forcetime;
|
|
|
|
|
dx[1]*=forcetime;
|
|
|
|
|
dx[2]*=forcetime;
|
|
|
|
|
@@ -2490,6 +2846,79 @@ static void softbody_restore_prev_step(Object *ob)
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
static void softbody_store_step(Object *ob)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there*/
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
int a;
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
VECCOPY(bp->prevvec,bp->vec);
|
|
|
|
|
VECCOPY(bp->prevpos,bp->pos);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* used by predictors and correctors */
|
|
|
|
|
static void softbody_store_state(Object *ob,float *ppos,float *pvel)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there*/
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
int a;
|
|
|
|
|
float *pp=ppos,*pv=pvel;
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
|
|
|
|
|
VECCOPY(pv, bp->vec);
|
|
|
|
|
pv+=3;
|
|
|
|
|
|
|
|
|
|
VECCOPY(pp, bp->pos);
|
|
|
|
|
pp+=3;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* used by predictors and correctors */
|
|
|
|
|
static void softbody_retrieve_state(Object *ob,float *ppos,float *pvel)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there*/
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
int a;
|
|
|
|
|
float *pp=ppos,*pv=pvel;
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
|
|
|
|
|
VECCOPY(bp->vec,pv);
|
|
|
|
|
pv+=3;
|
|
|
|
|
|
|
|
|
|
VECCOPY(bp->pos,pp);
|
|
|
|
|
pp+=3;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* used by predictors and correctors */
|
|
|
|
|
static void softbody_swap_state(Object *ob,float *ppos,float *pvel)
|
|
|
|
|
{
|
|
|
|
|
SoftBody *sb= ob->soft; /* is supposed to be there*/
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
int a;
|
|
|
|
|
float *pp=ppos,*pv=pvel;
|
|
|
|
|
float temp[3];
|
|
|
|
|
|
|
|
|
|
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
|
|
|
|
|
|
|
|
|
|
VECCOPY(temp, bp->vec);
|
|
|
|
|
VECCOPY(bp->vec,pv);
|
|
|
|
|
VECCOPY(pv,temp);
|
|
|
|
|
pv+=3;
|
|
|
|
|
|
|
|
|
|
VECCOPY(temp, bp->pos);
|
|
|
|
|
VECCOPY(bp->pos,pp);
|
|
|
|
|
VECCOPY(pp,temp);
|
|
|
|
|
pp+=3;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* care for bodypoints taken out of the 'ordinary' solver step
|
|
|
|
|
** because they are screwed to goal by bolts
|
|
|
|
|
** they just need to move along with the goal in time
|
|
|
|
|
@@ -3193,6 +3622,8 @@ SoftBody *sbNew(void)
|
|
|
|
|
|
|
|
|
|
sb->inspring= 0.5f;
|
|
|
|
|
sb->infrict= 0.5f;
|
|
|
|
|
/*todo backward file compat should copy infrict to inpush while reading old files*/
|
|
|
|
|
sb->inpush = 0.5f;
|
|
|
|
|
|
|
|
|
|
sb->interval= 10;
|
|
|
|
|
sb->sfra= G.scene->r.sfra;
|
|
|
|
|
@@ -3205,9 +3636,13 @@ SoftBody *sbNew(void)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
sb->minloops = 10;
|
|
|
|
|
sb->maxloops = 300;
|
|
|
|
|
|
|
|
|
|
sb->choke = 3;
|
|
|
|
|
sb_new_scratch(sb);
|
|
|
|
|
/*todo backward file compat should set sb->shearstiff = 1.0f while reading old files*/
|
|
|
|
|
sb->shearstiff = 1.0f;
|
|
|
|
|
sb->solverflags |= SBSO_OLDERR;
|
|
|
|
|
return sb;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
@@ -3256,7 +3691,7 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
HairKey *key= NULL;
|
|
|
|
|
BodyPoint *bp;
|
|
|
|
|
int a;
|
|
|
|
|
float dtime,ctime,forcetime,err;
|
|
|
|
|
float dtime,ctime,forcetime;
|
|
|
|
|
float hairmat[4][4];
|
|
|
|
|
|
|
|
|
|
/* This part only sets goals and springs, based on original mesh/curve/lattice data.
|
|
|
|
|
@@ -3446,8 +3881,10 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if (TRUE) { /* */
|
|
|
|
|
if (sb->solver_ID < 2) {
|
|
|
|
|
/* special case of 2nd order Runge-Kutta type AKA Heun */
|
|
|
|
|
int mid_flags=0;
|
|
|
|
|
float err = 0;
|
|
|
|
|
float forcetimemax = 1.0f;
|
|
|
|
|
float forcetimemin = 0.001f;
|
|
|
|
|
float timedone =0.0; /* how far did we get without violating error condition */
|
|
|
|
|
@@ -3459,7 +3896,8 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
|
|
|
|
|
|
|
|
|
|
if (sb->maxloops > 0) forcetimemin = 1.0f / sb->maxloops;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if(sb->solver_ID>0) mid_flags |= MID_PRESERVE;
|
|
|
|
|
|
|
|
|
|
//forcetime = dtime; /* hope for integrating in one step */
|
|
|
|
|
forcetime =forcetimemax; /* hope for integrating in one step */
|
|
|
|
|
@@ -3470,13 +3908,13 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
|
|
|
|
|
sb->scratch->flag &= ~SBF_DOFUZZY;
|
|
|
|
|
/* do predictive euler step */
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime);
|
|
|
|
|
softbody_apply_forces(ob, forcetime, 1, NULL);
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,0);
|
|
|
|
|
softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/* crop new slope values to do averaged slope step */
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime);
|
|
|
|
|
softbody_apply_forces(ob, forcetime, 2, &err);
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,0);
|
|
|
|
|
softbody_apply_forces(ob, forcetime, 2, &err,mid_flags);
|
|
|
|
|
|
|
|
|
|
softbody_apply_goalsnap(ob);
|
|
|
|
|
|
|
|
|
|
@@ -3517,6 +3955,7 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
sct=PIL_check_seconds_timer();
|
|
|
|
|
if (sct-sst > 0.5f) printf("%3.0f%% \r",100.0f*timedone);
|
|
|
|
|
}
|
|
|
|
|
/* ask for user break */
|
|
|
|
|
if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break;
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
@@ -3531,11 +3970,263 @@ void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
else if (sb->solver_ID == 2)
|
|
|
|
|
{
|
|
|
|
|
/* do semi "fake" implicit euler */
|
|
|
|
|
float *predict_vel=NULL,*predict_pos=NULL; /* for BDF style stepping */
|
|
|
|
|
NLContext *sb_nlc = NULL;
|
|
|
|
|
int npredict=0,predict_mem_size;
|
|
|
|
|
int nvar = 3*2*sb->totpoint;
|
|
|
|
|
int loops =0 ;
|
|
|
|
|
float forcetimemax = 1.0f; // this one needs 5 steps as a minimum
|
|
|
|
|
float forcetimemin = 0.001f;
|
|
|
|
|
float timedone =0.0; /* how far did we get without violating error condition */
|
|
|
|
|
/* loops = counter for emergency brake
|
|
|
|
|
* we don't want to lock up the system if physics fail
|
|
|
|
|
*/
|
|
|
|
|
SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */
|
|
|
|
|
if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
|
|
|
|
|
if (sb->maxloops > 0) forcetimemin = 1.0f / sb->maxloops;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
forcetime =forcetimemax; /* hope for integrating as fast as possible */
|
|
|
|
|
|
|
|
|
|
//allocate predictor buffers
|
|
|
|
|
npredict =1;
|
|
|
|
|
predict_mem_size =3*sizeof(float)*npredict*sb->totpoint;
|
|
|
|
|
predict_pos = MEM_mallocN(predict_mem_size,"SB_predict_pos");
|
|
|
|
|
predict_vel = MEM_mallocN(predict_mem_size,"SB_predict_vel");
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ){
|
|
|
|
|
float loss_of_accuracy=0;
|
|
|
|
|
// create new solver context for this loop
|
|
|
|
|
sb_nlc = (NLContext*)nlNewContext();
|
|
|
|
|
/* hum it smells like threading trouble */
|
|
|
|
|
nlSolverParameteri(NL_NB_VARIABLES, nvar);
|
|
|
|
|
nlSolverParameteri(NL_LEAST_SQUARES, NL_FALSE);
|
|
|
|
|
|
|
|
|
|
interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
|
|
|
|
|
softbody_store_step(ob); /* used for rolling back step guessing */
|
|
|
|
|
softbody_store_state(ob,predict_pos,predict_vel);
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,NLF_BUILD|NLF_SOLVE);
|
|
|
|
|
// go full step
|
|
|
|
|
softbody_apply_fake_implicit_forces(ob, forcetime, NULL,NULL,NULL);
|
|
|
|
|
|
|
|
|
|
// restore old situation and store 1rst solution to predictors
|
|
|
|
|
softbody_swap_state(ob,predict_pos,predict_vel);
|
|
|
|
|
// the following is to find out how good we were
|
|
|
|
|
// may be we can do smarter
|
|
|
|
|
// so now using the forces and jacobian we calculated before
|
|
|
|
|
// go only 1/2
|
|
|
|
|
softbody_apply_fake_implicit_forces(ob, forcetime/2.0f, NULL,NULL,NULL);
|
|
|
|
|
// explore situation here without redoing the jacobian
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,NLF_SOLVE);
|
|
|
|
|
// go next 1/2
|
|
|
|
|
softbody_apply_fake_implicit_forces(ob, forcetime/2.0f,&loss_of_accuracy,predict_pos,predict_vel );
|
|
|
|
|
// now we have "loss_of_accuracy"
|
|
|
|
|
|
|
|
|
|
softbody_apply_goalsnap(ob);
|
|
|
|
|
|
|
|
|
|
if (loss_of_accuracy > SoftHeunTol) { /* error needs to be scaled to some quantity */
|
|
|
|
|
|
|
|
|
|
if (forcetime > forcetimemin){
|
|
|
|
|
forcetime = MAX2(forcetime / 2.0f,forcetimemin);
|
|
|
|
|
softbody_restore_prev_step(ob);
|
|
|
|
|
//printf("down,");
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
timedone += forcetime;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
float newtime = forcetime * 1.5f; /* hope for 1.1 times better conditions in next step */
|
|
|
|
|
// all that 1/2 stepping was useless ... hum we know now ..
|
|
|
|
|
softbody_swap_state(ob,predict_pos,predict_vel);
|
|
|
|
|
if (0){//(sb->scratch->flag & SBF_DOFUZZY){
|
|
|
|
|
//if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */
|
|
|
|
|
newtime = forcetime;
|
|
|
|
|
//}
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
if (loss_of_accuracy > SoftHeunTol/1.1f) { /* stay with this stepsize unless err really small */
|
|
|
|
|
newtime = forcetime;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
timedone += forcetime;
|
|
|
|
|
newtime=MIN2(forcetimemax,MAX2(newtime,forcetimemin));
|
|
|
|
|
//if (newtime > forcetime) printf("up,");
|
|
|
|
|
if (forcetime > 0.0)
|
|
|
|
|
forcetime = MIN2(dtime - timedone,newtime);
|
|
|
|
|
else
|
|
|
|
|
forcetime = MAX2(dtime - timedone,newtime);
|
|
|
|
|
}
|
|
|
|
|
loops++;
|
|
|
|
|
// give away solver context within loop
|
|
|
|
|
if (sb_nlc)
|
|
|
|
|
{
|
|
|
|
|
if (sb_nlc != nlGetCurrent())printf("Aye NL context mismatch! in softbody.c !\n");
|
|
|
|
|
nlDeleteContext(sb_nlc);
|
|
|
|
|
sb_nlc = NULL;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if(sb->solverflags & SBSO_MONITOR ){
|
|
|
|
|
sct=PIL_check_seconds_timer();
|
|
|
|
|
if (sct-sst > 0.5f) printf("%3.0f%% \r",100.0f*timedone);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* ask for user break */
|
|
|
|
|
if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// give away buffers
|
|
|
|
|
if (predict_pos) MEM_freeN(predict_pos);
|
|
|
|
|
if (predict_vel) MEM_freeN(predict_vel);
|
|
|
|
|
|
|
|
|
|
if(sb->solverflags & SBSO_MONITOR ){
|
|
|
|
|
if (loops > HEUNWARNLIMIT) /* monitor high loop counts */
|
|
|
|
|
printf("\r needed %d steps/frame ",loops);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}/*SOLVER SELECT*/
|
|
|
|
|
else if (sb->solver_ID == 4)
|
|
|
|
|
{
|
|
|
|
|
/* do semi "fake" implicit euler */
|
|
|
|
|
NLContext *sb_nlc = NULL;
|
|
|
|
|
int nvar = 3*2*sb->totpoint;
|
|
|
|
|
int loops =0 ;
|
|
|
|
|
float forcetimemax = 1.0f; // this one needs 5 steps as a minimum
|
|
|
|
|
float forcetimemin = 0.001f;
|
|
|
|
|
float timedone =0.0; /* how far did we get without violating error condition */
|
|
|
|
|
/* loops = counter for emergency brake
|
|
|
|
|
* we don't want to lock up the system if physics fail
|
|
|
|
|
*/
|
|
|
|
|
SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */
|
|
|
|
|
if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
|
|
|
|
|
if (sb->maxloops > 0) forcetimemin = 1.0f / sb->maxloops;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
forcetime =forcetimemax; /* hope for integrating as fast as possible */
|
|
|
|
|
while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ){
|
|
|
|
|
float loss_of_accuracy=0;
|
|
|
|
|
// create new solver context for this loop
|
|
|
|
|
sb_nlc = (NLContext*)nlNewContext();
|
|
|
|
|
/* hum it smells like threading trouble */
|
|
|
|
|
nlSolverParameteri(NL_NB_VARIABLES, nvar);
|
|
|
|
|
nlSolverParameteri(NL_LEAST_SQUARES, NL_FALSE);
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
|
|
|
|
|
softbody_store_step(ob); /* used for rolling back step guessing */
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,NLF_BUILD|NLF_SOLVE);
|
|
|
|
|
softbody_apply_fake_implicit_forces(ob, forcetime,&loss_of_accuracy,NULL,NULL);
|
|
|
|
|
softbody_apply_goalsnap(ob);
|
|
|
|
|
|
|
|
|
|
if (loss_of_accuracy > SoftHeunTol) { /* error needs to be scaled to some quantity */
|
|
|
|
|
|
|
|
|
|
if (forcetime > forcetimemin){
|
|
|
|
|
forcetime = MAX2(forcetime / 2.0f,forcetimemin);
|
|
|
|
|
softbody_restore_prev_step(ob);
|
|
|
|
|
//printf("down,");
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
timedone += forcetime;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
float newtime = forcetime * 1.1f; /* hope for 1.1 times better conditions in next step */
|
|
|
|
|
if (0){//(sb->scratch->flag & SBF_DOFUZZY){
|
|
|
|
|
//if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */
|
|
|
|
|
newtime = forcetime;
|
|
|
|
|
//}
|
|
|
|
|
}
|
|
|
|
|
else {
|
|
|
|
|
if (loss_of_accuracy > SoftHeunTol/2.0f) { /* stay with this stepsize unless err really small */
|
|
|
|
|
newtime = forcetime;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
timedone += forcetime;
|
|
|
|
|
newtime=MIN2(forcetimemax,MAX2(newtime,forcetimemin));
|
|
|
|
|
//if (newtime > forcetime) printf("up,");
|
|
|
|
|
if (forcetime > 0.0)
|
|
|
|
|
forcetime = MIN2(dtime - timedone,newtime);
|
|
|
|
|
else
|
|
|
|
|
forcetime = MAX2(dtime - timedone,newtime);
|
|
|
|
|
}
|
|
|
|
|
loops++;
|
|
|
|
|
// give away solver context within loop
|
|
|
|
|
if (sb_nlc)
|
|
|
|
|
{
|
|
|
|
|
if (sb_nlc != nlGetCurrent())printf("Aye NL context mismatch! in softbody.c !\n");
|
|
|
|
|
nlDeleteContext(sb_nlc);
|
|
|
|
|
sb_nlc = NULL;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if(sb->solverflags & SBSO_MONITOR ){
|
|
|
|
|
sct=PIL_check_seconds_timer();
|
|
|
|
|
if (sct-sst > 0.5f) printf("%3.0f%% \r",100.0f*timedone);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* ask for user break */
|
|
|
|
|
if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
if(sb->solverflags & SBSO_MONITOR ){
|
|
|
|
|
if (loops > HEUNWARNLIMIT) /* monitor high loop counts */
|
|
|
|
|
printf("\r needed %d steps/frame ",loops);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}/*SOLVER SELECT*/
|
|
|
|
|
else if (sb->solver_ID == 3){
|
|
|
|
|
/* do "stupid" semi "fake" implicit euler */
|
|
|
|
|
NLContext *sb_nlc = NULL;
|
|
|
|
|
int nvar = 3*2*sb->totpoint;
|
|
|
|
|
int loops =0 ;
|
|
|
|
|
float forcetimemax = 1.0f; // this one needs 5 steps as a minimum
|
|
|
|
|
float timedone =0.0; /* how far did we get without violating error condition */
|
|
|
|
|
/* loops = counter for emergency brake
|
|
|
|
|
* we don't want to lock up the system if physics fail
|
|
|
|
|
*/
|
|
|
|
|
if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
forcetime =forcetimemax; /* hope for integrating as fast as possible */
|
|
|
|
|
|
|
|
|
|
while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ){
|
|
|
|
|
|
|
|
|
|
sb_nlc = (NLContext*)nlNewContext();
|
|
|
|
|
/* hum it smells like threading trouble */
|
|
|
|
|
nlSolverParameteri(NL_NB_VARIABLES, nvar);
|
|
|
|
|
nlSolverParameteri(NL_LEAST_SQUARES, NL_FALSE);
|
|
|
|
|
|
|
|
|
|
interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
|
|
|
|
|
softbody_calc_forces(ob, forcetime,timedone/dtime,NLF_BUILD|NLF_SOLVE);
|
|
|
|
|
softbody_apply_fake_implicit_forces(ob, forcetime, NULL,NULL,NULL);
|
|
|
|
|
softbody_apply_goalsnap(ob);
|
|
|
|
|
timedone += forcetime;
|
|
|
|
|
loops++;
|
|
|
|
|
if (sb_nlc)
|
|
|
|
|
{
|
|
|
|
|
if (sb_nlc != nlGetCurrent())printf("Aye NL context mismatch! in softbody.c !\n");
|
|
|
|
|
nlDeleteContext(sb_nlc);
|
|
|
|
|
sb_nlc = NULL;
|
|
|
|
|
}
|
|
|
|
|
/* ask for user break */
|
|
|
|
|
if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}/*SOLVER SELECT*/
|
|
|
|
|
else{
|
|
|
|
|
/* do brute force explicit euler */
|
|
|
|
|
/* removed but left this branch for better integrators / solvers (BM) */
|
|
|
|
|
/* yah! Nicholas Guttenberg (NichG) here is the place to plug in */
|
|
|
|
|
}
|
|
|
|
|
printf("softbody no valid solver ID!");
|
|
|
|
|
}/*SOLVER SELECT*/
|
|
|
|
|
|
|
|
|
|
if(sb->solverflags & SBSO_MONITOR ){
|
|
|
|
|
sct=PIL_check_seconds_timer();
|
|
|
|
|
if (sct-sst > 0.5f) printf(" solver time %f %s \r",sct-sst,ob->id.name);
|
|
|
|
|
|