diff options
Diffstat (limited to 'source/blender/blenkernel/intern/softbody.c')
| -rw-r--r-- | source/blender/blenkernel/intern/softbody.c | 2337 |
1 files changed, 1504 insertions, 833 deletions
diff --git a/source/blender/blenkernel/intern/softbody.c b/source/blender/blenkernel/intern/softbody.c index ac56d876df7..1ae547ba7b8 100644 --- a/source/blender/blenkernel/intern/softbody.c +++ b/source/blender/blenkernel/intern/softbody.c @@ -69,6 +69,8 @@ variables on the UI for now #include "BLI_blenlib.h" #include "BLI_arithb.h" #include "BLI_ghash.h" +#include "BLI_threads.h" + #include "BKE_curve.h" #include "BKE_effect.h" #include "BKE_global.h" @@ -84,6 +86,7 @@ variables on the UI for now //XXX #include "BIF_editdeform.h" //XXX #include "BIF_graphics.h" #include "PIL_time.h" +// #include "ONL_opennl.h" remove linking to ONL for now /* callbacks for errors and interrupts and some goo */ static int (*SB_localInterruptCallBack)(void) = NULL; @@ -94,7 +97,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; @@ -102,7 +105,7 @@ typedef struct BodySpring { typedef struct BodyFace { int v1, v2, v3 ,v4; - float ext_force[3]; /* edges colliding and sailing */ + float ext_force[3]; /* faces colliding */ short flag; } BodyFace; @@ -117,6 +120,25 @@ typedef struct SBScratch { float aabbmin[3],aabbmax[3]; }SBScratch; +typedef struct SB_thread_context{ + Object *ob; + float forcetime; + float timenow; + int ifirst; + int ilast; + ListBase *do_effector; + int do_deflector; + float fieldfactor; + float windfactor; + int nr; + int tot; +}SB_thread_context; + +#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 @@ -126,7 +148,9 @@ typedef struct SBScratch { #define BSF_INTERSECT 1 /* edge intersects collider face */ #define SBF_DOFUZZY 1 /* edge intersects collider face */ -#define BFF_INTERSECT 1 /* edge intersects collider face */ + +#define BFF_INTERSECT 1 /* collider edge intrudes face */ +#define BFF_CLOSEVERT 2 /* collider vertex repulses face */ float SoftHeunTol = 1.0f; /* humm .. this should be calculated from sb parameters and sizes */ @@ -212,7 +236,7 @@ typedef struct ccd_Mesh { -ccd_Mesh *ccd_mesh_make(Object *ob, DerivedMesh *dm) +static ccd_Mesh *ccd_mesh_make(Object *ob, DerivedMesh *dm) { ccd_Mesh *pccd_M = NULL; ccdf_minmax *mima =NULL; @@ -309,7 +333,7 @@ ccd_Mesh *ccd_mesh_make(Object *ob, DerivedMesh *dm) } return pccd_M; } -void ccd_mesh_update(Object *ob,ccd_Mesh *pccd_M, DerivedMesh *dm) +static void ccd_mesh_update(Object *ob,ccd_Mesh *pccd_M, DerivedMesh *dm) { ccdf_minmax *mima =NULL; MFace *mface=NULL; @@ -448,7 +472,7 @@ void ccd_mesh_update(Object *ob,ccd_Mesh *pccd_M, DerivedMesh *dm) return ; } -void ccd_mesh_free(ccd_Mesh *ccdm) +static void ccd_mesh_free(ccd_Mesh *ccdm) { if(ccdm && (ccdm->savety == CCD_SAVETY )){ /*make sure we're not nuking objects we don't know*/ MEM_freeN(ccdm->mface); @@ -460,7 +484,7 @@ void ccd_mesh_free(ccd_Mesh *ccdm) } } -void ccd_build_deflector_hache(Object *vertexowner,GHash *hash) +static void ccd_build_deflector_hache(Object *vertexowner,GHash *hash) { Base *base; Object *ob; @@ -512,7 +536,7 @@ void ccd_build_deflector_hache(Object *vertexowner,GHash *hash) } /* while (base) */ } -void ccd_update_deflector_hache(Object *vertexowner,GHash *hash) +static void ccd_update_deflector_hache(Object *vertexowner,GHash *hash) { Base *base; Object *ob; @@ -582,6 +606,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) { @@ -602,12 +627,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++; @@ -689,6 +714,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) { @@ -803,50 +829,6 @@ static void calculate_collision_balls(Object *ob) } -char set_octant_flags(float *ce, float *pos, float ball) -{ - float x,y,z; - char res = 0; - int a; - - for (a=0;a<7;a++){ - switch(a){ - case 0: x=pos[0]; y=pos[1]; z=pos[2]; break; - case 1: x=pos[0]+ball; y=pos[1]; z=pos[2]; break; - case 2: x=pos[0]-ball; y=pos[1]; z=pos[2]; break; - case 3: x=pos[0]; y=pos[1]+ball; z=pos[2]; break; - case 4: x=pos[0]; y=pos[1]-ball; z=pos[2]; break; - case 5: x=pos[0]; y=pos[1]; z=pos[2]+ball; break; - case 6: x=pos[0]; y=pos[1]; z=pos[2]-ball; break; - } - - x=pos[0]; y=pos[1]; z=pos[2]; - - if (x > ce[0]){ - if (y > ce[1]){ - if (z > ce[2]) res|= 1; - else res|= 2; - } - else{ - if (z > ce[2]) res|= 4; - else res|= 8; - } - } - - else{ - if (y > ce[1]){ - if (z > ce[2]) res|= 16; - else res|= 32; - } - else{ - if (z > ce[2]) res|= 64; - else res|= 128; - } - } - } - return res; -} - /* creates new softbody if didn't exist yet, makes new points and springs arrays */ static void renew_softbody(Object *ob, int totpoint, int totspring) { @@ -881,6 +863,8 @@ static void renew_softbody(Object *ob, int totpoint, int totspring) bp->nofsprings= 0; bp->springs= NULL; bp->choke = 0.0f; + bp->choke2 = 0.0f; + bp->frozen = 1.0f; bp->colball = 0.0f; bp->flag = 0; @@ -985,6 +969,7 @@ static void Vec3PlusStVec(float *v, float s, float *v1) } /* +++ dependancy information functions*/ + static int are_there_deflectors(unsigned int layer) { Base *base; @@ -1002,40 +987,11 @@ static int query_external_colliders(Object *me) { return(are_there_deflectors(me->lay)); } - -#if 0 -static int query_external_forces(Object *me) -{ -/* silly but true: we need to create effector cache to see if anything is in it */ - ListBase *ec = pdInitEffectors(me,NULL); - int result = 0; - if (ec){ - result = 1; - pdEndEffectors(ec); /* sorry ec, yes i'm an idiot, but i needed to know if you were there */ - } - return result; -} - -/* -any of that external objects may have an IPO or something alike .. -so unless we can ask them if they are moving we have to assume they do -*/ -static int query_external_time(Object *me) -{ - if (query_external_colliders(me)) return 1; - if (query_external_forces(me)) return 1; - return 0; -} -static int query_internal_time(Object *me) -{ - if (me->softflag & OB_SB_GOAL) return 1; - return 0; -} -#endif /* --- dependancy information functions*/ + /* +++ the aabb "force" section*/ -int sb_detect_aabb_collisionCached( float force[3], unsigned int par_layer,struct Object *vertexowner,float time) +static int sb_detect_aabb_collisionCached( float force[3], unsigned int par_layer,struct Object *vertexowner,float time) { Object *ob; SoftBody *sb=vertexowner->soft; @@ -1099,15 +1055,112 @@ int sb_detect_aabb_collisionCached( float force[3], unsigned int par_layer,struc /* +++ the face external section*/ +static int sb_detect_face_pointCached(float face_v1[3],float face_v2[3],float face_v3[3],float *damp, + float force[3], unsigned int par_layer,struct Object *vertexowner,float time) + { + Object *ob; + GHash *hash; + GHashIterator *ihash; + float nv1[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3],aabbmax[3]; + float facedist,outerfacethickness,tune = 10.f; + int a, deflected=0; + + aabbmin[0] = MIN3(face_v1[0],face_v2[0],face_v3[0]); + aabbmin[1] = MIN3(face_v1[1],face_v2[1],face_v3[1]); + aabbmin[2] = MIN3(face_v1[2],face_v2[2],face_v3[2]); + aabbmax[0] = MAX3(face_v1[0],face_v2[0],face_v3[0]); + aabbmax[1] = MAX3(face_v1[1],face_v2[1],face_v3[1]); + aabbmax[2] = MAX3(face_v1[2],face_v2[2],face_v3[2]); + + /* calculate face normal once again SIGH */ + VECSUB(edge1, face_v1, face_v2); + VECSUB(edge2, face_v3, face_v2); + Crossf(d_nvect, edge2, edge1); + Normalize(d_nvect); + -int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_v3[3],float *damp, + hash = vertexowner->soft->scratch->colliderhash; + ihash = BLI_ghashIterator_new(hash); + while (!BLI_ghashIterator_isDone(ihash) ) { + + ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); + ob = BLI_ghashIterator_getKey (ihash); + /* only with deflecting set */ + if(ob->pd && ob->pd->deflect) { + MVert *mvert= NULL; + MVert *mprevvert= NULL; + if(ccdm){ + mvert= ccdm->mvert; + a = ccdm->totvert; + mprevvert= ccdm->mprevvert; + outerfacethickness =ob->pd->pdef_sboft; + if ((aabbmax[0] < ccdm->bbmin[0]) || + (aabbmax[1] < ccdm->bbmin[1]) || + (aabbmax[2] < ccdm->bbmin[2]) || + (aabbmin[0] > ccdm->bbmax[0]) || + (aabbmin[1] > ccdm->bbmax[1]) || + (aabbmin[2] > ccdm->bbmax[2]) ) { + /* boxes dont intersect */ + BLI_ghashIterator_step(ihash); + continue; + } + + } + else{ + /*aye that should be cached*/ + printf("missing cache error \n"); + BLI_ghashIterator_step(ihash); + continue; + } + + + /* use mesh*/ + if (mvert) { + while(a){ + VECCOPY(nv1,mvert[a-1].co); + if(mprevvert){ + VecMulf(nv1,time); + Vec3PlusStVec(nv1,(1.0f-time),mprevvert[a-1].co); + } + /* origin to face_v2*/ + VECSUB(nv1, nv1, face_v2); + facedist = Inpf(nv1,d_nvect); + if (ABS(facedist)<outerfacethickness){ + if (point_in_tri_prism(nv1, face_v1,face_v2,face_v3) ){ + float df; + if (facedist > 0){ + df = (outerfacethickness-facedist)/outerfacethickness; + } + else { + df = (outerfacethickness+facedist)/outerfacethickness; + } + + *damp=df*tune*ob->pd->pdef_sbdamp; + + df = 0.01f*exp(- 100.0f*df); + Vec3PlusStVec(force,-df,d_nvect); + deflected = 3; + } + } + a--; + }/* while(a)*/ + } /* if (mvert) */ + } /* if(ob->pd && ob->pd->deflect) */ + BLI_ghashIterator_step(ihash); + } /* while () */ + BLI_ghashIterator_free(ihash); + return deflected; +} + + +static int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_v3[3],float *damp, float force[3], unsigned int par_layer,struct Object *vertexowner,float time) { Object *ob; GHash *hash; GHashIterator *ihash; float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3],aabbmax[3]; - float t; + float t,tune = 10.0f; int a, deflected=0; aabbmin[0] = MIN3(face_v1[0],face_v2[0],face_v3[0]); @@ -1206,8 +1259,8 @@ int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_ LineIntersectsTriangle(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) || LineIntersectsTriangle(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) || LineIntersectsTriangle(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL) ){ - Vec3PlusStVec(force,-1.0f,d_nvect); - *damp=ob->pd->pdef_sbdamp; + Vec3PlusStVec(force,-0.5f,d_nvect); + *damp=tune*ob->pd->pdef_sbdamp; deflected = 2; } if (mface->v4){ /* quad */ @@ -1217,11 +1270,12 @@ int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_ Crossf(d_nvect, edge2, edge1); Normalize(d_nvect); if ( - LineIntersectsTriangle(nv1, nv3, face_v1, face_v2, face_v3, &t, NULL) || + /* LineIntersectsTriangle(nv1, nv3, face_v1, face_v2, face_v3, &t, NULL) || + we did that edge allready */ LineIntersectsTriangle(nv3, nv4, face_v1, face_v2, face_v3, &t, NULL) || LineIntersectsTriangle(nv4, nv1, face_v1, face_v2, face_v3, &t, NULL) ){ - Vec3PlusStVec(force,-1.0f,d_nvect); - *damp=ob->pd->pdef_sbdamp; + Vec3PlusStVec(force,-0.5f,d_nvect); + *damp=tune*ob->pd->pdef_sbdamp; deflected = 2; } } @@ -1237,12 +1291,12 @@ int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_ -void scan_for_ext_face_forces(Object *ob,float timenow) +static void scan_for_ext_face_forces(Object *ob,float timenow) { SoftBody *sb = ob->soft; BodyFace *bf; int a; - float damp; + float damp=0.0f,choke=1.0f; float tune = -10.0f; float feedback[3]; @@ -1252,43 +1306,71 @@ void scan_for_ext_face_forces(Object *ob,float timenow) bf = sb->scratch->bodyface; for(a=0; a<sb->scratch->totface; a++, bf++) { bf->ext_force[0]=bf->ext_force[1]=bf->ext_force[2]=0.0f; +/*+++edges intruding*/ + bf->flag &= ~BFF_INTERSECT; feedback[0]=feedback[1]=feedback[2]=0.0f; - bf->flag &= ~BFF_INTERSECT; - if (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, &damp, feedback, ob->lay ,ob , timenow)){ - Vec3PlusStVec(bf->ext_force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v2].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); +// Vec3PlusStVec(bf->ext_force,tune,feedback); bf->flag |= BFF_INTERSECT; + choke = MIN2(MAX2(damp,choke),1.0f); } feedback[0]=feedback[1]=feedback[2]=0.0f; if ((bf->v4) && (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, &damp, feedback, ob->lay ,ob , timenow))){ - Vec3PlusStVec(bf->ext_force,tune,feedback); - bf->flag |= BFF_INTERSECT; + Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v4].force,tune,feedback); +// Vec3PlusStVec(bf->ext_force,tune,feedback); + bf->flag |= BFF_INTERSECT; + choke = MIN2(MAX2(damp,choke),1.0f); } - } +/*---edges intruding*/ + +/*+++ close vertices*/ + if (( bf->flag & BFF_INTERSECT)==0){ + bf->flag &= ~BFF_CLOSEVERT; + tune = -1.0f; + feedback[0]=feedback[1]=feedback[2]=0.0f; + if (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, + &damp, feedback, ob->lay ,ob , timenow)){ + Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v2].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); +// Vec3PlusStVec(bf->ext_force,tune,feedback); + bf->flag |= BFF_CLOSEVERT; + choke = MIN2(MAX2(damp,choke),1.0f); + } - - + feedback[0]=feedback[1]=feedback[2]=0.0f; + if ((bf->v4) && (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, + &damp, feedback, ob->lay ,ob , timenow))){ + Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); + Vec3PlusStVec(sb->bpoint[bf->v4].force,tune,feedback); +// Vec3PlusStVec(bf->ext_force,tune,feedback); + bf->flag |= BFF_CLOSEVERT; + choke = MIN2(MAX2(damp,choke),1.0f); + } + } +/*--- close vertices*/ + } bf = sb->scratch->bodyface; for(a=0; a<sb->scratch->totface; a++, bf++) { - if ( bf->flag & BFF_INTERSECT) + if (( bf->flag & BFF_INTERSECT) || ( bf->flag & BFF_CLOSEVERT)) { - VECADD(sb->bpoint[bf->v1].force,sb->bpoint[bf->v1].force,bf->ext_force); - VECADD(sb->bpoint[bf->v2].force,sb->bpoint[bf->v2].force,bf->ext_force); - VECADD(sb->bpoint[bf->v3].force,sb->bpoint[bf->v3].force,bf->ext_force); + sb->bpoint[bf->v1].choke2=MAX2(sb->bpoint[bf->v1].choke2,choke); + sb->bpoint[bf->v2].choke2=MAX2(sb->bpoint[bf->v2].choke2,choke); + sb->bpoint[bf->v3].choke2=MAX2(sb->bpoint[bf->v3].choke2,choke); if (bf->v4){ - VECADD(sb->bpoint[bf->v4].force,sb->bpoint[bf->v4].force,bf->ext_force); + sb->bpoint[bf->v2].choke2=MAX2(sb->bpoint[bf->v2].choke2,choke); } - - } - + } } - - - - } } @@ -1297,7 +1379,7 @@ void scan_for_ext_face_forces(Object *ob,float timenow) /* +++ the spring external section*/ -int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp, +static int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp, float force[3], unsigned int par_layer,struct Object *vertexowner,float time) { Object *ob; @@ -1447,17 +1529,16 @@ 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) + +static void _scan_for_ext_spring_forces(Object *ob,float timenow,int ifirst,int ilast, struct ListBase *do_effector) { SoftBody *sb = ob->soft; - ListBase *do_effector; int a; float damp; float feedback[3]; - do_effector= pdInitEffectors(ob,NULL); if (sb && sb->totspring){ - for(a=0; a<sb->totspring; a++) { + for(a=ifirst; a<ilast; a++) { BodySpring *bs = &sb->bspring[a]; bs->ext_force[0]=bs->ext_force[1]=bs->ext_force[2]=0.0f; feedback[0]=feedback[1]=feedback[2]=0.0f; @@ -1498,25 +1579,110 @@ void scan_for_ext_spring_forces(Object *ob,float timenow) } f = Normalize(vel); f = -0.0001f*f*f*sb->aeroedge; - /* todo add a nice angle dependant function */ - /* look up one at bergman scheafer */ + /* (todo) add a nice angle dependant function done for now BUT */ + /* still there could be some nice drag/lift function, but who needs it */ VECSUB(sp, sb->bpoint[bs->v1].pos , sb->bpoint[bs->v2].pos); Projf(pr,vel,sp); VECSUB(vel,vel,pr); Normalize(vel); - Vec3PlusStVec(bs->ext_force,f,vel); + if (ob->softflag & OB_SB_AERO_ANGLE){ + Normalize(sp); + Vec3PlusStVec(bs->ext_force,f*(1.0f-ABS(Inpf(vel,sp))),vel); + } + else{ + Vec3PlusStVec(bs->ext_force,f,vel); // to keep compatible with 2.45 release files + } } /* --- springs seeing wind */ } } } - if(do_effector) - pdEndEffectors(do_effector); } + + +static void scan_for_ext_spring_forces(Object *ob,float timenow) +{ + SoftBody *sb = ob->soft; + ListBase *do_effector= NULL; + do_effector= pdInitEffectors(ob,NULL); + if (sb){ + _scan_for_ext_spring_forces(ob,timenow,0,sb->totspring,do_effector); + } + if(do_effector) + pdEndEffectors(do_effector); +} + +static void *exec_scan_for_ext_spring_forces(void *data) +{ + SB_thread_context *pctx = (SB_thread_context*)data; + _scan_for_ext_spring_forces(pctx->ob,pctx->timenow,pctx->ifirst,pctx->ilast,pctx->do_effector); + return 0; +} + +static void sb_sfesf_threads_run(struct Object *ob, float timenow,int totsprings,int *ptr_to_break_func()) +{ + ListBase *do_effector = NULL; + ListBase threads; + SB_thread_context *sb_threads; + int i, totthread,left,dec; + int lowsprings =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ + + do_effector= pdInitEffectors(ob,NULL); + + /* figure the number of threads while preventing pretty pointless threading overhead */ + if(G.scene->r.mode & R_FIXED_THREADS) + totthread= G.scene->r.threads; + else + totthread= BLI_system_thread_count(); + /* what if we got zillions of CPUs running but less to spread*/ + while ((totsprings/totthread < lowsprings) && (totthread > 1)) { + totthread--; + } + + sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBSpringsThread"); + memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); + left = totsprings; + dec = totsprings/totthread +1; + for(i=0; i<totthread; i++) { + sb_threads[i].ob = ob; + sb_threads[i].forcetime = 0.0; // not used here + sb_threads[i].timenow = timenow; + sb_threads[i].ilast = left; + left = left - dec; + if (left >0){ + sb_threads[i].ifirst = left; + } + else + sb_threads[i].ifirst = 0; + sb_threads[i].do_effector = do_effector; + sb_threads[i].do_deflector = 0;// not used here + sb_threads[i].fieldfactor = 0.0f;// not used here + sb_threads[i].windfactor = 0.0f;// not used here + sb_threads[i].nr= i; + sb_threads[i].tot= totthread; + } + if(totthread > 1) { + BLI_init_threads(&threads, exec_scan_for_ext_spring_forces, totthread); + + for(i=0; i<totthread; i++) + BLI_insert_thread(&threads, &sb_threads[i]); + + BLI_end_threads(&threads); + } + else + exec_scan_for_ext_spring_forces(&sb_threads[0]); + /* clean up */ + MEM_freeN(sb_threads); + + if(do_effector) + pdEndEffectors(do_effector); +} + + /* --- the spring external section*/ -int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb,float*cc) +static int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb,float*cc) { float mindist,cp; int winner =1; @@ -1543,7 +1709,7 @@ int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb, -int sb_detect_vertex_collisionCached(float opco[3], float facenormal[3], float *damp, +static int sb_detect_vertex_collisionCached(float opco[3], float facenormal[3], float *damp, float force[3], unsigned int par_layer,struct Object *vertexowner, float time,float vel[3], float *intrusion) { @@ -1833,253 +1999,6 @@ int sb_detect_vertex_collisionCached(float opco[3], float facenormal[3], float * return deflected; } -/* not complete yet .. - try to find a pos resolving all inside collisions -*/ -#if 0 //mute it for now -int sb_detect_vertex_collisionCachedEx(float opco[3], float facenormal[3], float *damp, - float force[3], unsigned int par_layer,struct Object *vertexowner, - float time,float vel[3], float *intrusion) -{ - Object *ob; - GHash *hash; - GHashIterator *ihash; - float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3],d_nvect[3], dv1[3],ve[3],avel[3], - vv1[3], vv2[3], vv3[3], vv4[3], - facedist,n_mag,force_mag_norm,minx,miny,minz,maxx,maxy,maxz, - innerfacethickness,outerfacethickness, - closestinside, - ee = 5.0f, ff = 0.1f, fa; - int a, deflected=0, cavel=0; -/* init */ - *intrusion = 0.0f; - hash = vertexowner->soft->scratch->colliderhash; - ihash = BLI_ghashIterator_new(hash); -/* go */ - while (!BLI_ghashIterator_isDone(ihash) ) { - - ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); - ob = BLI_ghashIterator_getKey (ihash); - /* only with deflecting set */ - if(ob->pd && ob->pd->deflect) { - MFace *mface= NULL; - MVert *mvert= NULL; - MVert *mprevvert= NULL; - ccdf_minmax *mima= NULL; - - if(ccdm){ - mface= ccdm->mface; - mvert= ccdm->mvert; - mprevvert= ccdm->mprevvert; - mima= ccdm->mima; - a = ccdm->totface; - - minx =ccdm->bbmin[0]; - miny =ccdm->bbmin[1]; - minz =ccdm->bbmin[2]; - - maxx =ccdm->bbmax[0]; - maxy =ccdm->bbmax[1]; - maxz =ccdm->bbmax[2]; - - if ((opco[0] < minx) || - (opco[1] < miny) || - (opco[2] < minz) || - (opco[0] > maxx) || - (opco[1] > maxy) || - (opco[2] > maxz) ) { - /* outside the padded boundbox --> collision object is too far away */ - BLI_ghashIterator_step(ihash); - continue; - } - } - else{ - /*aye that should be cached*/ - printf("missing cache error \n"); - BLI_ghashIterator_step(ihash); - continue; - } - - /* do object level stuff */ - /* need to have user control for that since it depends on model scale */ - innerfacethickness =-ob->pd->pdef_sbift; - outerfacethickness =ob->pd->pdef_sboft; - closestinside = innerfacethickness; - fa = (ff*outerfacethickness-outerfacethickness); - fa *= fa; - fa = 1.0f/fa; - avel[0]=avel[1]=avel[2]=0.0f; - /* use mesh*/ - while (a--) { - if ( - (opco[0] < mima->minx) || - (opco[0] > mima->maxx) || - (opco[1] < mima->miny) || - (opco[1] > mima->maxy) || - (opco[2] < mima->minz) || - (opco[2] > mima->maxz) - ) { - mface++; - mima++; - continue; - } - - if (mvert){ - - VECCOPY(nv1,mvert[mface->v1].co); - VECCOPY(nv2,mvert[mface->v2].co); - VECCOPY(nv3,mvert[mface->v3].co); - if (mface->v4){ - VECCOPY(nv4,mvert[mface->v4].co); - } - - if (mprevvert){ - /* grab the average speed of the collider vertices - before we spoil nvX - humm could be done once a SB steps but then we' need to store that too - since the AABB reduced propabitlty to get here drasticallly - it might be a nice tradeof CPU <--> memory - */ - VECSUB(vv1,nv1,mprevvert[mface->v1].co); - VECSUB(vv2,nv2,mprevvert[mface->v2].co); - VECSUB(vv3,nv3,mprevvert[mface->v3].co); - if (mface->v4){ - VECSUB(vv4,nv4,mprevvert[mface->v4].co); - } - - VecMulf(nv1,time); - Vec3PlusStVec(nv1,(1.0f-time),mprevvert[mface->v1].co); - - VecMulf(nv2,time); - Vec3PlusStVec(nv2,(1.0f-time),mprevvert[mface->v2].co); - - VecMulf(nv3,time); - Vec3PlusStVec(nv3,(1.0f-time),mprevvert[mface->v3].co); - - if (mface->v4){ - VecMulf(nv4,time); - Vec3PlusStVec(nv4,(1.0f-time),mprevvert[mface->v4].co); - } - } - } - - /* switch origin to be nv2*/ - VECSUB(edge1, nv1, nv2); - VECSUB(edge2, nv3, nv2); - VECSUB(dv1,opco,nv2); /* abuse dv1 to have vertex in question at *origin* of triangle */ - - Crossf(d_nvect, edge2, edge1); - n_mag = Normalize(d_nvect); - facedist = Inpf(dv1,d_nvect); - - if ((facedist > closestinside) && (facedist < outerfacethickness)){ -// if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){ - if (point_in_tri_prism(opco, nv1, nv2, nv3) ){ - force_mag_norm =(float)exp(-ee*facedist); - if (facedist > outerfacethickness*ff) - force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); - *damp=ob->pd->pdef_sbdamp; - - if (facedist > 0.0f){ - *damp *= (1.0f - facedist/outerfacethickness); - Vec3PlusStVec(force,force_mag_norm,d_nvect); - if (deflected < 2){ - deflected = 1; - if ((mprevvert) && (*damp > 0.0f)){ - choose_winner(ve,opco,nv1,nv2,nv3,vv1,vv2,vv3); - VECADD(avel,avel,ve); - cavel ++; - } - } - - } - else{ - Vec3PlusStVec(force,force_mag_norm,d_nvect); - VECCOPY(facenormal,d_nvect); - if ((mprevvert) && (*damp > 0.0f)){ - choose_winner(avel,opco,nv1,nv2,nv3,vv1,vv2,vv3); - cavel = 1; - deflected = 2; - closestinside = facedist; - } - } - *intrusion = facedist; - } - } - if (mface->v4){ /* quad */ - /* switch origin to be nv4 */ - VECSUB(edge1, nv3, nv4); - VECSUB(edge2, nv1, nv4); - VECSUB(dv1,opco,nv4); /* abuse dv1 to have vertex in question at *origin* of triangle */ - - Crossf(d_nvect, edge2, edge1); - n_mag = Normalize(d_nvect); - facedist = Inpf(dv1,d_nvect); - - if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){ - if (point_in_tri_prism(opco, nv1, nv3, nv4) ){ - force_mag_norm =(float)exp(-ee*facedist); - if (facedist > outerfacethickness*ff) - force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); - Vec3PlusStVec(force,force_mag_norm,d_nvect); - *damp=ob->pd->pdef_sbdamp; - - if (facedist > 0.0f){ - *damp *= (1.0f - facedist/outerfacethickness); - Vec3PlusStVec(force,force_mag_norm,d_nvect); - if (deflected < 2){ - deflected = 1; - if ((mprevvert) && (*damp > 0.0f)){ - choose_winner(ve,opco,nv1,nv3,nv4,vv1,vv3,vv4); - VECADD(avel,avel,ve); - cavel ++; - } - } - - } - else{ - Vec3PlusStVec(force,force_mag_norm,d_nvect); - VECCOPY(facenormal,d_nvect); - if ((mprevvert) && (*damp > 0.0f)){ - choose_winner(avel,opco,nv1,nv3,nv4,vv1,vv3,vv4); - cavel = 1; - deflected = 2; - closestinside = facedist; - } - } - - - - *intrusion = facedist; - } - - } - } - mface++; - mima++; - }/* while a */ - } /* if(ob->pd && ob->pd->deflect) */ - BLI_ghashIterator_step(ihash); - } /* while () */ - BLI_ghashIterator_free(ihash); - if (cavel) VecMulf(avel,1.0f/(float)cavel); - VECCOPY(vel,avel); - - /* we did stay "outside" but have some close to contact forces - just to be complete fake a face normal - */ - if (deflected ==1){ - VECCOPY(facenormal,force); - Normalize(facenormal); - } - else{ - facenormal[0] = facenormal[1] = facenormal[2] = 0.0f; - } - return deflected; -} -#endif - - /* sandbox to plug in various deflection algos */ static int sb_deflect_face(Object *ob,float *actpos,float *facenormal,float *force,float *cf,float time,float *vel,float *intrusion) @@ -2092,102 +2011,191 @@ static int sb_deflect_face(Object *ob,float *actpos,float *facenormal,float *for return(deflected); } +/* hiding this for now .. but the jacobian may pop up on other tasks .. so i'd like to keep it +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 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 softbody_calc_forces(Object *ob, float forcetime, float timenow) +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) { -/* rule we never alter free variables :bp->vec bp->pos in here ! - * this will ruin adaptive stepsize AKA heun! (BM) - */ + 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 + } + } +} + + +/* since this is definitely the most CPU consuming task here .. try to spread it */ +/* core function _softbody_calc_forces_slice_in_a_thread */ +/* result is int to be able to flag user break */ +static int _softbody_calc_forces_slice_in_a_thread(Object *ob, float forcetime, float timenow,int ifirst,int ilast,int *ptr_to_break_func(),ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor) +{ + float iks; + int bb,do_selfcollision,do_springcollision,do_aero; + int number_of_points_here = ilast - ifirst; SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; - BodyPoint *bproot; - BodySpring *bs; - ListBase *do_effector; - float iks, ks, kd, gravity, actspringlen, forcefactor, sd[3]; - float fieldfactor = 1000.0f, windfactor = 250.0f; - float tune = sb->ballstiff; - int a, b, do_deflector,do_selfcollision,do_springcollision,do_aero; - - - gravity = sb->grav * sb_grav_force_scale(ob); - + /* intitialize */ + if (sb) { /* check conditions for various options */ - do_deflector= query_external_colliders(ob); - do_effector= pdInitEffectors(ob,NULL); + /* +++ could be done on object level to squeeze out the last bits of it */ do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); - - iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ - bproot= sb->bpoint; /* need this for proper spring addressing */ - - - - if (do_springcollision || do_aero) scan_for_ext_spring_forces(ob,timenow); - if (do_deflector) { - float defforce[3]; - do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow); + /* --- could be done on object level to squeeze out the last bits of it */ + } + else { + printf("Error expected a SB here \n"); + return (999); } - if (do_selfcollision ){ - float ce[3]; - VecMidf(ce,sb->scratch->aabbmax,sb->scratch->aabbmin); - for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { - bp->octantflag = set_octant_flags(ce,bp->pos,bp->colball); - } +/* debugerin */ + if (sb->totpoint < ifirst) { + printf("Aye 998"); + return (998); } +/* debugerin */ - for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { + + bp = &sb->bpoint[ifirst]; + for(bb=number_of_points_here; bb>0; bb--, bp++) { /* clear forces accumulator */ bp->force[0]= bp->force[1]= bp->force[2]= 0.0; - /* naive ball self collision */ /* needs to be done if goal snaps or not */ if(do_selfcollision){ int attached; BodyPoint *obp; + BodySpring *bs; int c,b; float velcenter[3],dvel[3],def[3]; float distance; float compare; + float bstune = sb->ballstiff; - for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) { - - if ((bp->octantflag & obp->octantflag) == 0) continue; - + for(c=sb->totpoint, obp= sb->bpoint; c>=ifirst+bb; c--, obp++) { compare = (obp->colball + bp->colball); VecSubf(def, bp->pos, obp->pos); - /* rather check the AABBoxes before ever calulating the real distance */ /* mathematically it is completly nuts, but performace is pretty much (3) times faster */ if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; - distance = Normalize(def); if (distance < compare ){ /* exclude body points attached with a spring */ attached = 0; for(b=obp->nofsprings;b>0;b--){ bs = sb->bspring + obp->springs[b-1]; - if (( sb->totpoint-a == bs->v2) || ( sb->totpoint-a == bs->v1)){ + if (( ilast-bb == bs->v2) || ( ilast-bb == bs->v1)){ attached=1; continue;} } if (!attached){ - float f = tune/(distance) + tune/(compare*compare)*distance - 2.0f*tune/compare ; + float f = bstune/(distance) + bstune/(compare*compare)*distance - 2.0f*bstune/compare ; VecMidf(velcenter, bp->vec, obp->vec); 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); + /* exploit force(a,b) == -force(b,a) part2/2 */ VecSubf(dvel,velcenter,obp->vec); VecMulf(dvel,sb->nodemass); Vec3PlusStVec(obp->force,sb->balldamp,dvel); Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def); - } } } @@ -2197,24 +2205,26 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow) if(bp->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */ 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); + float ks,kd; + 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]); + /* 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]); } else { bp->force[0]-= kd * (velgoal[0] - bp->vec[0]); @@ -2224,13 +2234,15 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow) } /* done goal stuff */ - /* gravitation */ + if (sb){ + float gravity = sb->grav * sb_grav_force_scale(ob); bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */ - + } /* particle field & vortex */ if(do_effector) { + float kd; float force[3]= {0.0f, 0.0f, 0.0f}; float speed[3]= {0.0f, 0.0f, 0.0f}; float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ @@ -2251,7 +2263,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow) } else { /* BP friction in media (not) moving*/ - kd= sb->mediafrict* sb_fric_force_scale(ob); + float kd = sb->mediafrict* sb_fric_force_scale(ob); /* assume it to be proportional to actual velocity */ bp->force[0]-= bp->vec[0]*kd; bp->force[1]-= bp->vec[1]*kd; @@ -2260,40 +2272,34 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow) } /* +++cached collision targets */ bp->choke = 0.0f; + bp->choke2 = 0.0f; bp->flag &= ~SBF_DOFUZZY; if(do_deflector) { float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion; - kd = 1.0f; + float kd = 1.0f; if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){ - if (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--> - we predict the solution beeing out of the collider - in heun step No1 and leave the heun step No2 adapt to it - so we kind of introduced a implicit solver for this case - */ - Vec3PlusStVec(bp->pos,-intrusion,facenormal); - - sb->scratch->flag |= SBF_DOFUZZY; - bp->flag |= SBF_DOFUZZY; - bp->choke = sb->choke*0.01f; - } - else{ + if (intrusion < 0.0f){ + sb->scratch->flag |= SBF_DOFUZZY; + bp->flag |= SBF_DOFUZZY; + bp->choke = sb->choke*0.01f; + } + VECSUB(cfforce,bp->vec,vel); Vec3PlusStVec(bp->force,-cf*50.0f,cfforce); - } - - Vec3PlusStVec(bp->force,kd,defforce); + + Vec3PlusStVec(bp->force,kd,defforce); } } /* ---cached collision targets */ /* +++springs */ + iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ if(ob->softflag & OB_SB_EDGES) { if (sb->bspring){ /* spring list exists at all ? */ + int b; + BodySpring *bs; for(b=bp->nofsprings;b>0;b--){ bs = sb->bspring + bp->springs[b-1]; if (do_springcollision || do_aero){ @@ -2302,67 +2308,532 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow) bp->choke = bs->cf; } - - if (( (sb->totpoint-a) == bs->v1) ){ - actspringlen= VecLenf( (bproot+bs->v2)->pos, bp->pos); - VecSubf(sd,(bproot+bs->v2)->pos, bp->pos); - 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) ){ - actspringlen= VecLenf( (bproot+bs->v1)->pos, bp->pos); - VecSubf(sd,bp->pos,(bproot+bs->v1)->pos); - 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,ilast-bb,bs,iks,forcetime,0); }/* loop springs */ }/* existing spring list */ }/*any edges*/ /* ---springs */ }/*omit on snap */ }/*loop all bp's*/ +return 0; /*done fine*/ +} + +static void *exec_softbody_calc_forces(void *data) +{ + SB_thread_context *pctx = (SB_thread_context*)data; + _softbody_calc_forces_slice_in_a_thread(pctx->ob,pctx->forcetime,pctx->timenow,pctx->ifirst,pctx->ilast,NULL,pctx->do_effector,pctx->do_deflector,pctx->fieldfactor,pctx->windfactor); + return 0; +} + +static void sb_cf_threads_run(struct Object *ob, float forcetime, float timenow,int totpoint,int *ptr_to_break_func(),struct ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor) +{ + ListBase threads; + SB_thread_context *sb_threads; + int i, totthread,left,dec; + int lowpoints =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ + + /* figure the number of threads while preventing pretty pointless threading overhead */ + if(G.scene->r.mode & R_FIXED_THREADS) + totthread= G.scene->r.threads; + else + totthread= BLI_system_thread_count(); + /* what if we got zillions of CPUs running but less to spread*/ + while ((totpoint/totthread < lowpoints) && (totthread > 1)) { + totthread--; + } + + /* printf("sb_cf_threads_run spawning %d threads \n",totthread); */ + + sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBThread"); + memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); + left = totpoint; + dec = totpoint/totthread +1; + for(i=0; i<totthread; i++) { + sb_threads[i].ob = ob; + sb_threads[i].forcetime = forcetime; + sb_threads[i].timenow = timenow; + sb_threads[i].ilast = left; + left = left - dec; + if (left >0){ + sb_threads[i].ifirst = left; + } + else + sb_threads[i].ifirst = 0; + sb_threads[i].do_effector = do_effector; + sb_threads[i].do_deflector = do_deflector; + sb_threads[i].fieldfactor = fieldfactor; + sb_threads[i].windfactor = windfactor; + sb_threads[i].nr= i; + sb_threads[i].tot= totthread; + } + + + if(totthread > 1) { + BLI_init_threads(&threads, exec_softbody_calc_forces, totthread); + + for(i=0; i<totthread; i++) + BLI_insert_thread(&threads, &sb_threads[i]); + + BLI_end_threads(&threads); + } + else + exec_softbody_calc_forces(&sb_threads[0]); + /* clean up */ + MEM_freeN(sb_threads); +} +static void softbody_calc_forcesEx(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) + */ + SoftBody *sb= ob->soft; /* is supposed to be there */ + BodyPoint *bproot; + ListBase *do_effector; + float iks, gravity; + float fieldfactor = 1000.0f, windfactor = 250.0f; + int do_deflector,do_selfcollision,do_springcollision,do_aero; + + gravity = sb->grav * sb_grav_force_scale(ob); + + /* check conditions for various options */ + do_deflector= query_external_colliders(ob); + do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); + do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); + do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); + + iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ + bproot= sb->bpoint; /* need this for proper spring addressing */ + + if (do_springcollision || do_aero) + sb_sfesf_threads_run(ob,timenow,sb->totspring,NULL); + + /* after spring scan because it uses Effoctors too */ + do_effector= pdInitEffectors(ob,NULL); + + if (do_deflector) { + float defforce[3]; + do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow); + } + + sb_cf_threads_run(ob,forcetime,timenow,sb->totpoint,NULL,do_effector,do_deflector,fieldfactor,windfactor); /* finally add forces caused by face collision */ - if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow); - /* cleanup */ + if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow); + + /* finish matrix and solve */ if(do_effector) pdEndEffectors(do_effector); } -static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err) + +static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags) +{ + /* redirection to the new threaded Version */ + if (!(G.rt & 0x10)){ // 16 + softbody_calc_forcesEx(ob, forcetime, timenow, nl_flags); + return; + } + else{ + /* so the following will die */ + /* |||||||||||||||||||||||||| */ + /* VVVVVVVVVVVVVVVVVVVVVVVVVV */ + /*backward compatibility note: + fixing bug [17428] which forces adaptive step size to tiny steps + in some situations + .. keeping G.rt==17 0x11 option for old files 'needing' the bug*/ + + /* rule we never alter free variables :bp->vec bp->pos in here ! + * this will ruin adaptive stepsize AKA heun! (BM) + */ + SoftBody *sb= ob->soft; /* is supposed to be there */ + BodyPoint *bp; + BodyPoint *bproot; + BodySpring *bs; + ListBase *do_effector; + 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; + + + /* jacobian + NLboolean success; + + if(nl_flags){ + nlBegin(NL_SYSTEM); + nlBegin(NL_MATRIX); + } + */ + + + gravity = sb->grav * sb_grav_force_scale(ob); + + /* check conditions for various options */ + do_deflector= query_external_colliders(ob); + do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); + do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); + do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); + + iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ + bproot= sb->bpoint; /* need this for proper spring addressing */ + + if (do_springcollision || do_aero) scan_for_ext_spring_forces(ob,timenow); + /* after spring scan because it uses Effoctors too */ + do_effector= pdInitEffectors(ob,NULL); + + if (do_deflector) { + float defforce[3]; + do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow); + } + + 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 */ + /* jacobioan + nlMatrixAdd(op+ia,ia,-forcetime); + nlMatrixAdd(op+ia+1,ia+1,-forcetime); + nlMatrixAdd(op+ia+2,ia+2,-forcetime); + + 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 */ + if(do_selfcollision){ + int attached; + BodyPoint *obp; + int c,b; + float velcenter[3],dvel[3],def[3]; + float distance; + float compare; + + for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) { + + //if ((bp->octantflag & obp->octantflag) == 0) continue; + + compare = (obp->colball + bp->colball); + VecSubf(def, bp->pos, obp->pos); + + /* rather check the AABBoxes before ever calulating the real distance */ + /* mathematically it is completly nuts, but performace is pretty much (3) times faster */ + if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; + + distance = Normalize(def); + if (distance < compare ){ + /* exclude body points attached with a spring */ + attached = 0; + for(b=obp->nofsprings;b>0;b--){ + bs = sb->bspring + obp->springs[b-1]; + if (( sb->totpoint-a == bs->v2) || ( sb->totpoint-a == bs->v1)){ + attached=1; + continue;} + } + if (!attached){ + float f = tune/(distance) + tune/(compare*compare)*distance - 2.0f*tune/compare ; + + VecMidf(velcenter, bp->vec, obp->vec); + VecSubf(dvel,velcenter,bp->vec); + VecMulf(dvel,sb->nodemass); + + 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); + + Vec3PlusStVec(obp->force,sb->balldamp,dvel); + Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def); + + + } + } + } + } + /* naive ball self collision done */ + + if(bp->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */ + float auxvect[3]; + float velgoal[3]; + + /* do goal stuff */ + if(ob->softflag & OB_SB_GOAL) { + /* true elastic goal */ + 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]); + + 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 * (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]); + bp->force[1]-= kd * (velgoal[1] - bp->vec[1]); + bp->force[2]-= kd * (velgoal[2] - bp->vec[2]); + } + } + /* done goal stuff */ + + + /* 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 */ + if(do_effector) { + float force[3]= {0.0f, 0.0f, 0.0f}; + float speed[3]= {0.0f, 0.0f, 0.0f}; + float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ + + pdDoEffectors(do_effector, bp->pos, force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED); + + /* apply forcefield*/ + VecMulf(force,fieldfactor* eval_sb_fric_force_scale); + VECADD(bp->force, bp->force, force); + + /* BP friction in moving media */ + kd= sb->mediafrict* eval_sb_fric_force_scale; + bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale); + bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale); + bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale); + /* now we'll have nice centrifugal effect for vortex */ + + } + else { + /* BP friction in media (not) moving*/ + kd= sb->mediafrict* sb_fric_force_scale(ob); + /* assume it to be proportional to actual velocity */ + bp->force[0]-= bp->vec[0]*kd; + 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); + /* 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; + bp->choke2 = 0.0f; + bp->flag &= ~SBF_DOFUZZY; + if(do_deflector) { + float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion; + kd = 1.0f; + + if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){ + if ((!nl_flags)&&(intrusion < 0.0f)){ + if(G.rt & 0x01){ // 17 we did check for bit 0x10 before + /*fixing bug [17428] this forces adaptive step size to tiny steps + in some situations .. keeping G.rt==17 option for old files 'needing' the bug + */ + /*bjornmose: uugh.. what an evil hack + violation of the 'don't touch bp->pos in here' rule + but works nice, like this--> + we predict the solution beeing out of the collider + in heun step No1 and leave the heun step No2 adapt to it + so we kind of introduced a implicit solver for this case + */ + Vec3PlusStVec(bp->pos,-intrusion,facenormal); + } + else{ + + VECSUB(cfforce,bp->vec,vel); + Vec3PlusStVec(bp->force,-cf*50.0f,cfforce); + } + + + sb->scratch->flag |= SBF_DOFUZZY; + bp->flag |= SBF_DOFUZZY; + bp->choke = sb->choke*0.01f; + } + else{ + 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); + + } + + } + + } + /* ---cached collision targets */ + + /* +++springs */ + if(ob->softflag & OB_SB_EDGES) { + if (sb->bspring){ /* spring list exists at all ? */ + for(b=bp->nofsprings;b>0;b--){ + bs = sb->bspring + bp->springs[b-1]; + if (do_springcollision || do_aero){ + VecAddf(bp->force,bp->force,bs->ext_force); + if (bs->flag & BSF_INTERSECT) + bp->choke = bs->cf; + + } + // sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags) + // rather remove nl_falgs from code .. will make things a lot cleaner + sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,0); + }/* loop springs */ + }/* existing spring list */ + }/*any edges*/ + /* ---springs */ + }/*omit on snap */ + }/*loop all bp's*/ + + + /* finally add forces caused by face collision */ + if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow); + + /* finish matrix and solve */ +#if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be usefull .. so keep that BM + 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 == 32) && (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 ==32) + 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 */ +#endif + if(do_effector) pdEndEffectors(do_effector); + } +} + +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 */ @@ -2370,7 +2841,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float * SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; float dx[3],dv[3],aabbmin[3],aabbmax[3],cm[3]={0.0f,0.0f,0.0f}; - float timeovermass; + float timeovermass/*,freezeloc=0.00001f,freezeforce=0.00000000001f*/; float maxerrpos= 0.0f,maxerrvel = 0.0f; int a,fuzzy=0; @@ -2385,14 +2856,14 @@ 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){ + /* this makes t~ = t */ + 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 */ /* the euler step for velocity then becomes */ /* v(t + dt) = v(t) + a(t) * dt */ - bp->force[0]*= timeovermass; /* individual mass of node here */ - bp->force[1]*= timeovermass; - bp->force[2]*= timeovermass; + VecMulf(bp->force,timeovermass);/* individual mass of node here */ /* some nasty if's to have heun in here too */ VECCOPY(dv,bp->force); @@ -2413,14 +2884,24 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float * } else {VECADD(bp->vec, bp->vec, bp->force);} + /* this makes t~ = t+dt */ + if(!(mid_flags & MID_PRESERVE)) VECCOPY(dx,bp->vec); + /* 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; + /* x(t + dt) = x(t) + v(t~) * dt */ + VecMulf(dx,forcetime); + + /* the freezer coming sooner or later */ + /* + if ((Inpf(dx,dx)<freezeloc )&&(Inpf(bp->force,bp->force)<freezeforce )){ + bp->frozen /=2; + } + else{ + bp->frozen =MIN2(bp->frozen*1.05f,1.0f); + } + VecMulf(dx,bp->frozen); + */ /* again some nasty if's to have heun in here too */ if (mode ==1){ VECCOPY(bp->prevpos,bp->pos); @@ -2430,7 +2911,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float * if (mode ==2){ bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]); bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]); - bp->pos[2] = bp->prevpos[2] + 0.5f* ( dx[2] + bp->prevdx[2]); + bp->pos[2] = bp->prevpos[2] + 0.5f * ( dx[2] + bp->prevdx[2]); maxerrpos = MAX2(maxerrpos,ABS(dx[0] - bp->prevdx[0])); maxerrpos = MAX2(maxerrpos,ABS(dx[1] - bp->prevdx[1])); maxerrpos = MAX2(maxerrpos,ABS(dx[2] - bp->prevdx[2])); @@ -2439,10 +2920,11 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float * 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->choke2 > 0.0f){ + VecMulf(bp->vec,(1.0f - bp->choke2)); + } if (bp->choke > 0.0f){ - bp->vec[0] = bp->vec[0]*(1.0f - bp->choke); - bp->vec[1] = bp->vec[1]*(1.0f - bp->choke); - bp->vec[2] = bp->vec[2]*(1.0f - bp->choke); + VecMulf(bp->vec,(1.0f - bp->choke)); } } @@ -2489,6 +2971,81 @@ static void softbody_restore_prev_step(Object *ob) } } +#if 0 +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; + } +} +#endif + + /* 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 @@ -2509,6 +3066,30 @@ static void softbody_apply_goalsnap(Object *ob) } } + +static void apply_spring_memory(Object *ob) +{ + SoftBody *sb = ob->soft; + BodySpring *bs; + BodyPoint *bp1, *bp2; + int a; + float b,l,r; + + b = sb->plastic; + if (sb && sb->totspring){ + for(a=0; a<sb->totspring; a++) { + bs = &sb->bspring[a]; + bp1 =&sb->bpoint[bs->v1]; + bp2 =&sb->bpoint[bs->v2]; + l = VecLenf(bp1->pos,bp2->pos); + r = bs->len/l; + if (( r > 1.05f) || (r < 0.95)){ + bs->len = ((100.0f - b) * bs->len + b*l)/100.0f; + } + } + } +} + /* expects full initialized softbody */ static void interpolate_exciter(Object *ob, int timescale, int time) { @@ -2597,9 +3178,9 @@ static void springs_from_mesh(Object *ob) } /* recalculate spring length for meshes here */ - /* special hidden feature! shrink to fit */ - if (G.rt > 500){ - scale = (G.rt - 500) / 100.0f; + /* public version shrink to fit */ + if (sb->springpreload != 0 ){ + scale = sb->springpreload / 100.0f; } for(a=0; a<sb->totspring; a++) { BodySpring *bs = &sb->bspring[a]; @@ -2712,7 +3293,7 @@ static void mesh_faces_to_scratch(Object *ob) helper function to get proper spring length when object is rescaled */ -float globallen(float *v1,float *v2,Object *ob) +static float globallen(float *v1,float *v2,Object *ob) { float p1[3],p2[3]; VECCOPY(p1,v1); @@ -3002,9 +3583,12 @@ static void particles_to_softbody(Object *ob) /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ renew_softbody(ob, totpoint, totedge); - psys->particles->bpi = 0; - for(a=1, pa=psys->particles+1; a<psys->totpart; a++, pa++) - pa->bpi = (pa-1)->bpi + pa->totkey; + /* find first BodyPoint index for each particle */ + if(psys->totpart > 0) { + psys->particles->bpi = 0; + for(a=1, pa=psys->particles+1; a<psys->totpart; a++, pa++) + pa->bpi = (pa-1)->bpi + (pa-1)->totkey; + } /* we always make body points */ sb= ob->soft; @@ -3066,97 +3650,68 @@ static void softbody_to_object(Object *ob, float (*vertexCos)[3], int numVerts, } } -void softbody_clear_cache(Object *ob, float framenr) -{ - SoftBody *sb = ob->soft; - ModifierData *md = ob->modifiers.first; - int stack_index = -1; - int a; - - if(sb==NULL) return; - - if(sb->particles) - stack_index = modifiers_indexInObject(ob,(ModifierData*)psys_get_modifier(ob,sb->particles)); - else { - for(a=0; md; md=md->next, a++) { - if(md->type == eModifierType_Softbody) { - stack_index = a; - break; - } - } - } - - BKE_ptcache_id_clear((ID *)ob, PTCACHE_CLEAR_ALL, framenr, stack_index); -} -static void softbody_write_cache(Object *ob, float framenr) +void sbWriteCache(Object *ob, int framenr) { - FILE *fp = NULL; - SoftBody *sb = ob->soft; + SoftBody *sb= ob->soft; BodyPoint *bp; - ModifierData *md = ob->modifiers.first; - int stack_index = -1; + PTCacheID pid; + PTCacheFile *pf; int a; - if(sb->totpoint == 0) return; - - if(sb->particles) - stack_index = modifiers_indexInObject(ob,(ModifierData*)psys_get_modifier(ob,sb->particles)); - else { - for(a=0; md; md=md->next, a++) { - if(md->type == eModifierType_Softbody) { - stack_index = a; - break; - } - } - } + if(sb->totpoint == 0) + return; - fp = BKE_ptcache_id_fopen((ID *)ob, 'w', framenr, stack_index); - if(!fp) return; + BKE_ptcache_id_from_softbody(&pid, ob, sb); + pf= BKE_ptcache_file_open(&pid, PTCACHE_FILE_WRITE, framenr); + if(!pf) + return; + + for(a=0, bp=sb->bpoint; a<sb->totpoint; a++, bp++) + BKE_ptcache_file_write_floats(pf, bp->pos, 3); for(a=0, bp=sb->bpoint; a<sb->totpoint; a++, bp++) - fwrite(&bp->pos, sizeof(float), 3, fp); - - fclose(fp); + BKE_ptcache_file_write_floats(pf, bp->vec, 3); + + BKE_ptcache_file_close(pf); } + static int softbody_read_cache(Object *ob, float framenr) { - FILE *fp = NULL; - SoftBody *sb = ob->soft; + SoftBody *sb= ob->soft; BodyPoint *bp; - ModifierData *md = ob->modifiers.first; - int stack_index = -1; - int a, ret = 1; + PTCacheID pid; + PTCacheFile *pf; + int a; - if(sb->totpoint == 0) return 0; + if(sb->totpoint == 0) + return 0; + + BKE_ptcache_id_from_softbody(&pid, ob, sb); + pf= BKE_ptcache_file_open(&pid, PTCACHE_FILE_READ, framenr); + if(!pf) + return 0; - if(sb->particles) - stack_index = modifiers_indexInObject(ob,(ModifierData*)psys_get_modifier(ob,sb->particles)); - else { - for(a=0; md; md=md->next, a++) { - if(md->type == eModifierType_Softbody) { - stack_index = a; - break; - } + for(a=0, bp=sb->bpoint; a<sb->totpoint; a++, bp++) { + if(!BKE_ptcache_file_read_floats(pf, bp->pos, 3)) { + BKE_ptcache_file_close(pf); + return 0; } } - fp = BKE_ptcache_id_fopen((ID *)ob, 'r', framenr, stack_index); - if(!fp) - ret = 0; - else { - for(a=0, bp=sb->bpoint; a<sb->totpoint; a++, bp++) - if(fread(&bp->pos, sizeof(float), 3, fp) != 3) { - ret = 0; - break; - } - - fclose(fp); + for(a=0, bp=sb->bpoint; a<sb->totpoint; a++, bp++) { + if(!BKE_ptcache_file_read_floats(pf, bp->vec, 3)) { + BKE_ptcache_file_close(pf); + return 0; + } } - return ret; + BKE_ptcache_file_close(pf); + + return 1; } + /* +++ ************ maintaining scratch *************** */ -void sb_new_scratch(SoftBody *sb) +static void sb_new_scratch(SoftBody *sb) { if (!sb) return; sb->scratch = MEM_callocN(sizeof(SBScratch), "SBScratch"); @@ -3192,6 +3747,8 @@ SoftBody *sbNew(void) sb->inspring= 0.5f; sb->infrict= 0.5f; + /*todo backward file compat should copy inspring to inpush while reading old files*/ + sb->inpush = 0.5f; sb->interval= 10; sb->sfra= G.scene->r.sfra; @@ -3204,9 +3761,16 @@ 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; + + sb->pointcache = BKE_ptcache_add(); + return sb; } @@ -3214,14 +3778,19 @@ SoftBody *sbNew(void) void sbFree(SoftBody *sb) { free_softbody_intern(sb); + BKE_ptcache_free(sb->pointcache); MEM_freeN(sb); } +void sbFreeSimulation(SoftBody *sb) +{ + free_softbody_intern(sb); +} /* makes totally fresh start situation */ void sbObjectToSoftbody(Object *ob) { - ob->softflag |= OB_SB_REDO; + //ob->softflag |= OB_SB_REDO; free_softbody_intern(ob->soft); } @@ -3245,306 +3814,408 @@ void sbSetInterruptCallBack(int (*f)(void)) SB_localInterruptCallBack = f; } - -/* simulates one step. framenr is in frames */ -void sbObjectStep(Object *ob, float framenr, float (*vertexCos)[3], int numVerts) +static void softbody_update_positions(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) { - ParticleSystemModifierData *psmd=0; - ParticleData *pa=0; - SoftBody *sb; + ParticleSystemModifierData *psmd= NULL; + ParticleData *pa= NULL; HairKey *key= NULL; BodyPoint *bp; + float hairmat[4][4]; int a; - float dtime,ctime,forcetime,err; + + /* update the vertex locations */ + if(sb->particles && sb->particles->totpart>0) { + psmd= psys_get_modifier(ob,sb->particles); + + pa= sb->particles->particles; + key= pa->hair; + + psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); + } + + for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { + /* store where goals are now */ + VECCOPY(bp->origS, bp->origE); + /* copy the position of the goals at desired end time */ + if(sb->particles) { + if(key == pa->hair + pa->totkey) { + pa++; + key = pa->hair; + + psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); + } + VECCOPY(bp->origE, key->co); + Mat4MulVecfl(hairmat,bp->origE); + + key++; + } + else{ + VECCOPY(bp->origE, vertexCos[a]); + /* vertexCos came from local world, go global */ + Mat4MulVecfl(ob->obmat, bp->origE); + } + /* just to be save give bp->origT a defined value + will be calulated in interpolate_exciter()*/ + VECCOPY(bp->origT, bp->origE); + } +} + +static void softbody_reset(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) +{ + ParticleSystemModifierData *psmd= NULL; + HairKey *key= NULL; + ParticleData *pa= NULL; + BodyPoint *bp; float hairmat[4][4]; + int a; - /* This part only sets goals and springs, based on original mesh/curve/lattice data. - Copying coordinates happens in next chunk by setting softbody flag OB_SB_RESET */ - /* remake softbody if: */ - if( (ob->softflag & OB_SB_REDO) || /* signal after weightpainting */ - (ob->soft==NULL) || /* just to be nice we allow full init */ - (ob->soft->bpoint==NULL) || /* after reading new file, or acceptable as signal to refresh */ - (numVerts!=ob->soft->totpoint) || /* should never happen, just to be safe */ - ((ob->softflag & OB_SB_EDGES) && !ob->soft->bspring && object_has_edges(ob))) /* happens when in UI edges was set */ - { - if(ob->soft && ob->soft->bpoint) /* don't clear on file load */ - softbody_clear_cache(ob, framenr); + if(sb->particles && sb->particles->totpart>0) { + psmd= psys_get_modifier(ob, sb->particles); + pa= sb->particles->particles; + key= pa->hair; - if(ob->soft->particles){ - particles_to_softbody(ob); + psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); + } + + for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { + if(sb->particles) { + if(key == pa->hair + pa->totkey) { + pa++; + key = pa->hair; + + psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); + } + VECCOPY(bp->pos, key->co); + Mat4MulVecfl(hairmat, bp->pos); + key++; + } + else { + VECCOPY(bp->pos, vertexCos[a]); + Mat4MulVecfl(ob->obmat, bp->pos); /* yep, sofbody is global coords*/ } - else switch(ob->type) { + VECCOPY(bp->origS, bp->pos); + VECCOPY(bp->origE, bp->pos); + VECCOPY(bp->origT, bp->pos); + bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f; + + /* the bp->prev*'s are for rolling back from a canceled try to propagate in time + adaptive step size algo in a nutshell: + 1. set sheduled time step to new dtime + 2. try to advance the sheduled time step, beeing optimistic execute it + 3. check for success + 3.a we 're fine continue, may be we can increase sheduled time again ?? if so, do so! + 3.b we did exceed error limit --> roll back, shorten the sheduled time and try again at 2. + 4. check if we did reach dtime + 4.a nope we need to do some more at 2. + 4.b yup we're done + */ + + VECCOPY(bp->prevpos, bp->pos); + VECCOPY(bp->prevvec, bp->vec); + VECCOPY(bp->prevdx, bp->vec); + VECCOPY(bp->prevdv, bp->vec); + } + + /* make a nice clean scratch struc */ + free_scratch(sb); /* clear if any */ + sb_new_scratch(sb); /* make a new */ + sb->scratch->needstobuildcollider=1; + + if((sb->particles)==0) { + /* copy some info to scratch */ + switch(ob->type) { case OB_MESH: - mesh_to_softbody(ob); + if (ob->softflag & OB_SB_FACECOLL) mesh_faces_to_scratch(ob); break; case OB_LATTICE: - lattice_to_softbody(ob); break; case OB_CURVE: case OB_SURF: - curve_surf_to_softbody(ob); break; default: - renew_softbody(ob, numVerts, 0); break; } + } +} + +static void softbody_step(Object *ob, SoftBody *sb, float dtime) +{ + /* the simulator */ + float forcetime; + double sct,sst=PIL_check_seconds_timer(); + ccd_update_deflector_hache(ob,sb->scratch->colliderhash); - /* still need to update to correct vertex locations, happens on next step */ - ob->softflag |= OB_SB_RESET; - ob->softflag &= ~OB_SB_REDO; + if(sb->scratch->needstobuildcollider){ + if (query_external_colliders(ob)){ + ccd_build_deflector_hache(ob,sb->scratch->colliderhash); + } + sb->scratch->needstobuildcollider=0; } - sb= ob->soft; + 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 */ + /* loops = counter for emergency brake + * we don't want to lock up the system if physics fail + */ + int loops =0 ; + + 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; - /* still no points? go away */ - if(sb->totpoint==0) return; + 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 */ + while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ) + { + /* set goals in time */ + interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime))); + + sb->scratch->flag &= ~SBF_DOFUZZY; + /* do predictive euler step */ + softbody_calc_forces(ob, forcetime,timedone/dtime,0); - if(sb->particles){ - psmd=psys_get_modifier(ob,sb->particles); - pa=sb->particles->particles; - } + softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags); - /* checking time: */ + /* crop new slope values to do averaged slope step */ + softbody_calc_forces(ob, forcetime,timedone/dtime,0); - ctime= bsystem_time(ob, framenr, 0.0); + softbody_apply_forces(ob, forcetime, 2, &err,mid_flags); + softbody_apply_goalsnap(ob); + + if (err > 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 (sb->scratch->flag & SBF_DOFUZZY){ + //if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */ + newtime = forcetime; + //} + } + else { + if (err > 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++; + 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 (ob->softflag&OB_SB_RESET) { - dtime = 0.0; - } else { - dtime= ctime - sb->ctime; + } + /* move snapped to final position */ + interpolate_exciter(ob, 2, 2); + softbody_apply_goalsnap(ob); + + // if(G.f & G_DEBUG){ + if(sb->solverflags & SBSO_MONITOR ){ + if (loops > HEUNWARNLIMIT) /* monitor high loop counts */ + printf("\r needed %d steps/frame",loops); + } + } + else if (sb->solver_ID == 2) + {/* do semi "fake" implicit euler */ + //removed + }/*SOLVER SELECT*/ + else if (sb->solver_ID == 4) + { + /* do semi "fake" implicit euler */ + }/*SOLVER SELECT*/ + else if (sb->solver_ID == 3){ + /* do "stupid" semi "fake" implicit euler */ + //removed - if(softbody_read_cache(ob, framenr)) { - if(sb->particles==0) - softbody_to_object(ob, vertexCos, numVerts, sb->local); - sb->ctime = ctime; - return; + }/*SOLVER SELECT*/ + else{ + printf("softbody no valid solver ID!"); + }/*SOLVER SELECT*/ + if(sb->plastic){ apply_spring_memory(ob);} + + if(sb->solverflags & SBSO_MONITOR ){ + sct=PIL_check_seconds_timer(); + if ((sct-sst > 0.5f) || (G.f & G_DEBUG)) printf(" solver time %f sec %s \n",sct-sst,ob->id.name); } +} - /* the simulator */ +/* simulates one step. framenr is in frames */ +void sbObjectStep(Object *ob, float cfra, float (*vertexCos)[3], int numVerts) +{ + ParticleSystemModifierData *psmd=0; + ParticleData *pa=0; + SoftBody *sb= ob->soft; + PointCache *cache; + PTCacheID pid; + float dtime, timescale; + int framedelta, framenr, startframe, endframe; - /* update the vertex locations */ - if (dtime!=0.0) { - if(sb->particles) { - pa=sb->particles->particles; - key = pa->hair; + cache= sb->pointcache; - psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); - } + framenr= (int)cfra; + framedelta= framenr - cache->simframe; - for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { - /* store where goals are now */ - VECCOPY(bp->origS, bp->origE); - /* copy the position of the goals at desired end time */ - if(sb->particles) { - if(key == pa->hair + pa->totkey) { - pa++; - key = pa->hair; + BKE_ptcache_id_from_softbody(&pid, ob, sb); + BKE_ptcache_id_time(&pid, framenr, &startframe, &endframe, ×cale); - psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); - } - VECCOPY(bp->origE, key->co); - Mat4MulVecfl(hairmat,bp->origE); + /* check for changes in mesh, should only happen in case the mesh + * structure changes during an animation */ + if(sb->bpoint && numVerts != sb->totpoint) { + cache->flag &= ~PTCACHE_SIMULATION_VALID; + cache->simframe= 0; - key++; - } - else{ - VECCOPY(bp->origE, vertexCos[a]); - /* vertexCos came from local world, go global */ - Mat4MulVecfl(ob->obmat, bp->origE); - } - /* just to be save give bp->origT a defined value - will be calulated in interpolate_exciter()*/ - VECCOPY(bp->origT, bp->origE); - } + return; } - if((ob->softflag&OB_SB_RESET) || /* got a reset signal */ - (dtime<0.0) || /* back in time */ - (dtime>=9.9*G.scene->r.framelen) /* too far forward in time --> goals won't be accurate enough */ - ) - { - if(sb->particles) { - pa=sb->particles->particles; - key = pa->hair; + /* clamp frame ranges */ + if(framenr < startframe) { + cache->flag &= ~PTCACHE_SIMULATION_VALID; + cache->simframe= 0; - psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); - } + return; + } + else if(framenr > endframe) { + framenr = endframe; + } - for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { - if(sb->particles) { - if(key == pa->hair + pa->totkey) { - pa++; - key = pa->hair; + /* verify if we need to create the softbody data */ + if(sb->bpoint == NULL || + ((ob->softflag & OB_SB_EDGES) && !ob->soft->bspring && object_has_edges(ob))) { - psys_mat_hair_to_global(ob, psmd->dm, sb->particles->part->from, pa, hairmat); - } - VECCOPY(bp->pos, key->co); - Mat4MulVecfl(hairmat, bp->pos); - key++; - } - else { - VECCOPY(bp->pos, vertexCos[a]); - Mat4MulVecfl(ob->obmat, bp->pos); /* yep, sofbody is global coords*/ - } - VECCOPY(bp->origS, bp->pos); - VECCOPY(bp->origE, bp->pos); - VECCOPY(bp->origT, bp->pos); - bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f; - - /* the bp->prev*'s are for rolling back from a canceled try to propagate in time - adaptive step size algo in a nutshell: - 1. set sheduled time step to new dtime - 2. try to advance the sheduled time step, beeing optimistic execute it - 3. check for success - 3.a we 're fine continue, may be we can increase sheduled time again ?? if so, do so! - 3.b we did exceed error limit --> roll back, shorten the sheduled time and try again at 2. - 4. check if we did reach dtime - 4.a nope we need to do some more at 2. - 4.b yup we're done - */ - - VECCOPY(bp->prevpos, bp->pos); - VECCOPY(bp->prevvec, bp->vec); - VECCOPY(bp->prevdx, bp->vec); - VECCOPY(bp->prevdv, bp->vec); + if(sb->particles){ + particles_to_softbody(ob); } - /* make a nice clean scratch struc */ - free_scratch(sb); /* clear if any */ - sb_new_scratch(sb); /* make a new */ - sb->scratch->needstobuildcollider=1; - - if((sb->particles)==0) { - /* copy some info to scratch */ + else { switch(ob->type) { - case OB_MESH: - if (ob->softflag & OB_SB_FACECOLL) mesh_faces_to_scratch(ob); - break; - case OB_LATTICE: - break; - case OB_CURVE: - case OB_SURF: - break; - default: - break; + case OB_MESH: + mesh_to_softbody(ob); + break; + case OB_LATTICE: + lattice_to_softbody(ob); + break; + case OB_CURVE: + case OB_SURF: + curve_surf_to_softbody(ob); + break; + default: + renew_softbody(ob, numVerts, 0); + break; } } - ob->softflag &= ~OB_SB_RESET; + softbody_update_positions(ob, sb, vertexCos, numVerts); + softbody_reset(ob, sb, vertexCos, numVerts); } - else if(dtime>0.0) { - double sct,sst=PIL_check_seconds_timer(); - ccd_update_deflector_hache(ob,sb->scratch->colliderhash); + /* continue physics special case */ + if(BKE_ptcache_get_continue_physics()) { + cache->flag &= ~PTCACHE_SIMULATION_VALID; + cache->simframe= 0; - if(sb->scratch->needstobuildcollider){ - if (query_external_colliders(ob)){ - ccd_build_deflector_hache(ob,sb->scratch->colliderhash); - } - sb->scratch->needstobuildcollider=0; - } + /* do simulation */ + dtime = timescale; + softbody_update_positions(ob, sb, vertexCos, numVerts); + softbody_step(ob, sb, dtime); - if (TRUE) { /* */ - /* special case of 2nd order Runge-Kutta type AKA Heun */ - float forcetimemax = 1.0f; - 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 - */ - int loops =0 ; - 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 = dtime; /* hope for integrating in one step */ - forcetime =forcetimemax; /* hope for integrating in one step */ - while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ) - { - /* set goals in time */ - interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime))); - - sb->scratch->flag &= ~SBF_DOFUZZY; - /* do predictive euler step */ - softbody_calc_forces(ob, forcetime,timedone/dtime); - softbody_apply_forces(ob, forcetime, 1, NULL); + if(sb->particles==0) + softbody_to_object(ob, vertexCos, numVerts, 0); + return; + } - /* crop new slope values to do averaged slope step */ - softbody_calc_forces(ob, forcetime,timedone/dtime); - softbody_apply_forces(ob, forcetime, 2, &err); + /* still no points? go away */ + if(sb->totpoint==0) return; - softbody_apply_goalsnap(ob); - - if (err > 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 (sb->scratch->flag & SBF_DOFUZZY){ - //if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */ - newtime = forcetime; - //} - } - else { - if (err > 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++; - if(sb->solverflags & SBSO_MONITOR ){ - sct=PIL_check_seconds_timer(); - if (sct-sst > 0.5f) printf("%3.0f%% \r",100.0f*timedone); - } - if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break; + if(sb->particles){ + psmd= psys_get_modifier(ob, sb->particles); + pa= sb->particles->particles; + } - } - /* move snapped to final position */ - interpolate_exciter(ob, 2, 2); - softbody_apply_goalsnap(ob); - - // if(G.f & G_DEBUG){ - if(sb->solverflags & SBSO_MONITOR ){ - if (loops > HEUNWARNLIMIT) /* monitor high loop counts */ - printf("\r needed %d steps/frame ",loops); - } - - } - 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 */ - } - 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); - } + /* try to read from cache */ + if(softbody_read_cache(ob, framenr)) { + if(sb->particles==0) + softbody_to_object(ob, vertexCos, numVerts, sb->local); + + cache->flag |= PTCACHE_SIMULATION_VALID; + cache->simframe= framenr; + + return; + } + else if(ob->id.lib || (cache->flag & PTCACHE_BAKED)) { + /* if baked and nothing in cache, do nothing */ + if(cache->flag & PTCACHE_SIMULATION_VALID) { + cache->flag &= ~PTCACHE_SIMULATION_VALID; + cache->simframe= 0; + } + + return; } - if(sb->particles==0) - softbody_to_object(ob, vertexCos, numVerts, 0); - sb->ctime= ctime; + if(framenr == startframe) { + /* first frame, no simulation to do, just set the positions */ + softbody_update_positions(ob, sb, vertexCos, numVerts); + + cache->flag |= PTCACHE_SIMULATION_VALID; + cache->simframe= framenr; + + /* don't write cache on first frame, but on second frame write + * cache for frame 1 and 2 */ + } + else if(framedelta == 1) { + /* if on second frame, write cache for first frame */ + if(framenr == startframe+1) + sbWriteCache(ob, startframe); + + softbody_update_positions(ob, sb, vertexCos, numVerts); + + /* do simulation */ + cache->flag |= PTCACHE_SIMULATION_VALID; + cache->simframe= framenr; + + /* checking time: */ + dtime = framedelta*timescale; - softbody_write_cache(ob, framenr); + softbody_step(ob, sb, dtime); + + if(sb->particles==0) + softbody_to_object(ob, vertexCos, numVerts, 0); + + sbWriteCache(ob, framenr); + } + else { + /* time step backwards or too large forward - do nothing */ + if(cache->flag & PTCACHE_SIMULATION_VALID) { + cache->flag &= ~PTCACHE_SIMULATION_VALID; + cache->simframe= 0; + } + } } |