[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [13177] trunk/blender/source/blender: quite a load is still hidden behind the define #ifdef _work_on_sb_solver
Jens Ole Wund (bjornmose)
bjornmose at gmx.net
Wed Jan 9 01:25:53 CET 2008
Revision: 13177
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=13177
Author: bjornmose
Date: 2008-01-09 01:25:51 +0100 (Wed, 09 Jan 2008)
Log Message:
-----------
quite a load is still hidden behind the define #ifdef _work_on_sb_solver
a glance to view is the "STU PID semi implicit euler"
most of the work to implement a semi implicit euler was done ..
now i am dealing with the tradeoffs between 'calculation time' which is quite expensive .. inverting a 0(n*n) sparse matrix /* once agian thanks to brecht for his work on making sparse matrices LU decomposition and evaluating inverses that easy*/ putting it into and cropping works pretty nice .. argh off topic again
...
while i spent a little time on reading recent papers i found :
1. control on springs needs to be split in pushing and pulling /* fabric pushes easy but pulls hard */
2. diagonals on 4-gons (in the current SB model) can be seen as shear .. thus need a contol to modify .. this commit wil add it
3. 2 nd order springs /*aka rigidity */ can focus on bending .. thus renaming 'em
i have no idea how i would provide backward compatiblity, but the spots i marked in code :)
Modified Paths:
--------------
trunk/blender/source/blender/blenkernel/BKE_softbody.h
trunk/blender/source/blender/blenkernel/SConscript
trunk/blender/source/blender/blenkernel/intern/softbody.c
trunk/blender/source/blender/makesdna/DNA_object_force.h
trunk/blender/source/blender/src/buttons_object.c
Modified: trunk/blender/source/blender/blenkernel/BKE_softbody.h
===================================================================
--- trunk/blender/source/blender/blenkernel/BKE_softbody.h 2008-01-09 00:05:26 UTC (rev 13176)
+++ trunk/blender/source/blender/blenkernel/BKE_softbody.h 2008-01-09 00:25:51 UTC (rev 13177)
@@ -41,6 +41,7 @@
float origS[3], origE[3], origT[3], pos[3], vec[3], force[3];
float goal;
float prevpos[3], prevvec[3], prevdx[3], prevdv[3]; /* used for Heun integration */
+ float impdv[3],impdx[3];
int nofsprings; int *springs;
float choke;
float colball;
Modified: trunk/blender/source/blender/blenkernel/SConscript
===================================================================
--- trunk/blender/source/blender/blenkernel/SConscript 2008-01-09 00:05:26 UTC (rev 13176)
+++ trunk/blender/source/blender/blenkernel/SConscript 2008-01-09 00:25:51 UTC (rev 13177)
@@ -8,6 +8,7 @@
incs += ' ../imbuf ../avi #/intern/elbeem/extern ../nodes'
incs += ' #/intern/iksolver/extern ../blenloader ../quicktime'
incs += ' #/intern/bmfont'
+incs += ' #/intern/opennl/extern'
incs += ' ' + env['BF_OPENGL_INC']
incs += ' ' + env['BF_ZLIB_INC']
Modified: trunk/blender/source/blender/blenkernel/intern/softbody.c
===================================================================
--- trunk/blender/source/blender/blenkernel/intern/softbody.c 2008-01-09 00:05:26 UTC (rev 13176)
+++ trunk/blender/source/blender/blenkernel/intern/softbody.c 2008-01-09 00:25:51 UTC (rev 13177)
@@ -87,6 +87,7 @@
#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 @@
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 @@
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 @@
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 @@
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 @@
{
int counter = 0;
BodySpring *bs_new;
+ stiffness *=stiffness;
add_2nd_order_roller(ob,stiffness,&counter,0); /* counting */
if (counter) {
@@ -1450,6 +1458,7 @@
}
+
void scan_for_ext_spring_forces(Object *ob,float timenow)
{
SoftBody *sb = ob->soft;
@@ -2095,9 +2104,111 @@
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 @@
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 @@
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,11 +2260,34 @@
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 */
if(do_selfcollision){
@@ -2156,7 +2300,7 @@
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 @@
VecSubf(dvel,velcenter,bp->vec);
VecMulf(dvel,sb->nodemass);
+ Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
Vec3PlusStVec(bp->force,sb->balldamp,dvel);
- Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
+
+ 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 @@
Vec3PlusStVec(obp->force,sb->balldamp,dvel);
Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def);
+
}
}
}
@@ -2201,23 +2371,39 @@
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);
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list