source/blender/blenkernel/intern/softbody.c

   1 /*  softbody.c
   2  *
   3  * $Id:
   4  *
   5  * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version. The Blender
  11  * Foundation also sells licenses for use in proprietary software under
  12  * the Blender License.  See http://www.blender.org/BL/ for information
  13  * about this.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License
  21  * along with this program; if not, write to the Free Software Foundation,
  22  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  23  *
  24  * The Original Code is Copyright (C) Blender Foundation
  25  * All rights reserved.
  26  *
  27  * The Original Code is: all of this file.
  28  *
  29  * Contributor(s): none yet.
  30  *
  31  * ***** END GPL/BL DUAL LICENSE BLOCK *****
  32  */
  33
  34 /*
  35 ******
  36 variables on the UI for now
  37
  38         float mediafrict;  friction to env
  39         float nodemass;   softbody mass of *vertex*
  40         float grav;        softbody amount of gravitaion to apply
  41
  42         float goalspring;  softbody goal springs
  43         float goalfrict;   softbody goal springs friction
  44         float mingoal;     quick limits for goal
  45         float maxgoal;
  46
  47         float inspring;   softbody inner springs
  48         float infrict;     softbody inner springs friction
  49
  50 *****
  51 */
  52
  53
  54 #include <math.h>
  55 #include <stdlib.h>
  56 #include <string.h>
  57
  58 #include "MEM_guardedalloc.h"
  59
  60 /* types */
  61 #include "DNA_curve_types.h"
  62 #include "DNA_object_types.h"
  63 #include "DNA_mesh_types.h"
  64 #include "DNA_meshdata_types.h"
  65 #include "DNA_lattice_types.h"
  66 #include "DNA_scene_types.h"
  67
  68 #include "BLI_blenlib.h"
  69 #include "BLI_arithb.h"
  70
  71 #include "BKE_displist.h"
  72 #include "BKE_effect.h"
  73 #include "BKE_global.h"
  74 #include "BKE_object.h"
  75 #include "BKE_softbody.h"
  76 #include "BKE_utildefines.h"
  77
  78 #include  "BIF_editdeform.h"
  79
  80 extern bDeformGroup *get_named_vertexgroup(Object *ob, char *name);
  81 extern int  get_defgroup_num (Object *ob, bDeformGroup        *dg);
  82
  83
  84 /* ********** soft body engine ******* */
  85 #define SOFTGOALSNAP  0.999f
  86 // if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
  87 // removes *unnecessary* stiffnes from ODE system
  88 #define HEUNWARNLIMIT 1 // 50 would be fine i think for detecting severe *stiff* stuff
  89
  90
  91 float SoftHeunTol = 1.0f; // humm .. this should be calculated from sb parameters and sizes
  92
  93 /* local prototypes */
  94 static void free_softbody_intern(SoftBody *sb);
  95
  96
  97 /*+++ frame based timing +++*/
  98
  99 //physical unit of force is [kg * m / sec^2]
 100
 101 static float sb_grav_force_scale(Object *ob)
 102 // since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm
 103 // put it to a function here, so we can add user options later without touching simulation code
 104 {
 105         return (0.001f);
 106 }
 107
 108 static float sb_fric_force_scale(Object *ob)
 109 // rescaling unit of drag [1 / sec] to somehow reasonable
 110 // put it to a function here, so we can add user options later without touching simulation code
 111 {
 112         return (0.01f);
 113 }
 114
 115 static float sb_time_scale(Object *ob)
 116 // defining the frames to *real* time relation
 117 {
 118         SoftBody *sb= ob->soft; // is supposed to be there
 119         if (sb){
 120                 return(sb->physics_speed); //hrms .. this could be IPO as well :)
 121                 // estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)]
 122                 // 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions]  has period 65 frames
 123         // theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM)
 124         }
 125         return (1.0f);
 126         /*
 127         this would be frames/sec independant timing assuming 25 fps is default
 128         but does not work very well with NLA
 129                 return (25.0f/G.scene->r.frs_sec)
 130         */
 131 }
 132 /*--- frame based timing ---*/
 133
 134
 135 static int count_mesh_quads(Mesh *me)
 136 {
 137         int a,result = 0;
 138         MFace *mface= me->mface;
 139
 140         if(mface) {
 141                 for(a=me->totface; a>0; a--, mface++) {
 142                         if(mface->v4) result++;
 143                 }
 144         }
 145         return result;
 146 }
 147
 148 static void add_mesh_quad_diag_springs(Object *ob)
 149 {
 150         Mesh *me= ob->data;
 151         MFace *mface= me->mface;
 152         BodyPoint *bp;
 153         BodySpring *bs, *bs_new;
 154         int a ;
 155
 156         if (ob->soft){
 157                 int nofquads;
 158
 159                 nofquads = count_mesh_quads(me);
 160                 if (nofquads) {
 161                         /* resize spring-array to hold additional quad springs */
 162                         bs_new= MEM_callocN( (ob->soft->totspring + nofquads *2 )*sizeof(BodySpring), "bodyspring");
 163                         memcpy(bs_new,ob->soft->bspring,(ob->soft->totspring )*sizeof(BodySpring));
 164
 165                         if(ob->soft->bspring)
 166                                 MEM_freeN(ob->soft->bspring); /* do this before reassigning the pointer  or have a 1st class memory leak */
 167                         ob->soft->bspring = bs_new;
 168
 169                         /* fill the tail */
 170                         a = 0;
 171                         bs = bs_new+ob->soft->totspring;
 172                         bp= ob->soft->bpoint;
 173                         if(mface ) {
 174                                 for(a=me->totface; a>0; a--, mface++) {
 175                                         if(mface->v4) {
 176                                                 bs->v1= mface->v1;
 177                                                 bs->v2= mface->v3;
 178                                                 bs->strength= 1.0;
 179                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
 180                                                 bs++;
 181                                                 bs->v1= mface->v2;
 182                                                 bs->v2= mface->v4;
 183                                                 bs->strength= 1.0;
 184                                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
 185                                                 bs++;
 186
 187                                         }
 188                                 }
 189                         }
 190
 191             /* now we can announce new springs */
 192                         ob->soft->totspring += nofquads *2;
 193                 }
 194         }
 195 }
 196
 197
 198 static void add_bp_springlist(BodyPoint *bp,int springID)
 199 {
 200         int *newlist;
 201
 202         if (bp->springs == NULL) {
 203                 bp->springs = MEM_callocN( sizeof(int), "bpsprings");
 204                 bp->springs[0] = springID;
 205                 bp->nofsprings = 1;
 206         }
 207         else {
 208                 bp->nofsprings++;
 209                 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings");
 210                 memcpy(newlist,bp->springs,(bp->nofsprings-1)* sizeof(int));
 211                 MEM_freeN(bp->springs);
 212                 bp->springs = newlist;
 213                 bp->springs[bp->nofsprings-1] = springID;
 214         }
 215 }
 216
 217 /* do this once when sb is build
 218 it is O(N^2) so scanning for springs every iteration is too expensive
 219 */
 220 static void build_bps_springlist(Object *ob)
 221 {
 222         SoftBody *sb= ob->soft; // is supposed to be there
 223         BodyPoint *bp;
 224         BodySpring *bs;
 225         int a,b;
 226
 227         if (sb==NULL) return; // paranoya check
 228
 229         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 230                 /* scan for attached inner springs */
 231                 for(b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) {
 232                         if (( (sb->totpoint-a) == bs->v1) ){
 233                                 add_bp_springlist(bp,sb->totspring -b);
 234                         }
 235                         if (( (sb->totpoint-a) == bs->v2) ){
 236                                 add_bp_springlist(bp,sb->totspring -b);
 237                         }
 238                 }//for springs
 239                 // if (bp->nofsprings) printf(" node %d has %d spring links\n",a,bp->nofsprings);
 240         }//for bp
 241 }
 242
 243
 244 /* creates new softbody if didn't exist yet, makes new points and springs arrays */
 245 /* called in mesh_to_softbody */
 246 static void renew_softbody(Object *ob, int totpoint, int totspring)
 247 {
 248         SoftBody *sb;
 249
 250         if(ob->soft==NULL) ob->soft= sbNew();
 251         else free_softbody_intern(ob->soft);
 252         sb= ob->soft;
 253
 254         if(totpoint) {
 255                 sb->totpoint= totpoint;
 256                 sb->totspring= totspring;
 257
 258                 sb->bpoint= MEM_mallocN( totpoint*sizeof(BodyPoint), "bodypoint");
 259                 if(totspring)
 260                         sb->bspring= MEM_mallocN( totspring*sizeof(BodySpring), "bodyspring");
 261         }
 262 }
 263
 264 /* only frees internal data */
 265 static void free_softbody_intern(SoftBody *sb)
 266 {
 267         if(sb) {
 268                 int a;
 269                 BodyPoint *bp;
 270
 271                 if(sb->bpoint){
 272                         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 273                                 /* free spring list */
 274                                 if (bp->springs != NULL) {
 275                                         MEM_freeN(bp->springs);
 276                                 }
 277                         }
 278                         MEM_freeN(sb->bpoint);
 279                 }
 280
 281                 if(sb->bspring) MEM_freeN(sb->bspring);
 282
 283                 sb->totpoint= sb->totspring= 0;
 284                 sb->bpoint= NULL;
 285                 sb->bspring= NULL;
 286         }
 287 }
 288
 289
 290 /* ************ dynamics ********** */
 291
 292 /* aye this belongs to arith.c */
 293 static void Vec3PlusStVec(float *v, float s, float *v1)
 294 {
 295         v[0] += s*v1[0];
 296         v[1] += s*v1[1];
 297         v[2] += s*v1[2];
 298 }
 299
 300 static int sb_deflect_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *force,float *cf ,float *bounce)
 301 {
 302         int deflected;
 303         float s_actpos[3], s_futurepos[3];
 304         VECCOPY(s_actpos,actpos);
 305         if(futurepos)
 306         VECCOPY(s_futurepos,futurepos);
 307         if (bounce) *bounce *= 1.5f;
 308
 309
 310         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
 311                                         facenormal, cf, force , 1,
 312                                         G.scene->r.cfra, ob->lay, ob);
 313         return(deflected);
 314
 315 }
 316
 317 /* for future use (BM)
 318 static int sb_deflect_edge_face(Object *ob,float *actpos, float *futurepos,float *collisionpos, float *facenormal,float *slip ,float *bounce)
 319 {
 320         int deflected;
 321         float dummy[3],s_actpos[3], s_futurepos[3];
 322         SoftBody *sb= ob->soft; // is supposed to be there
 323         VECCOPY(s_actpos,actpos);
 324         VECCOPY(s_futurepos,futurepos);
 325         if (slip)   *slip   *= 0.98f;
 326         if (bounce) *bounce *= 1.5f;
 327
 328
 329         deflected= SoftBodyDetectCollision(s_actpos, s_futurepos, collisionpos,
 330                                         facenormal, dummy, dummy , 2,
 331                                         G.scene->r.cfra, ob->lay, ob);
 332         return(deflected);
 333
 334 }
 335 */
 336 // some functions removed here .. to help HOS on next merge (BM)
 337
 338 #define USES_FIELD              1
 339 #define USES_DEFLECT    2
 340 static int is_there_deflection(unsigned int layer)
 341 {
 342         Base *base;
 343         int retval= 0;
 344
 345         for(base = G.scene->base.first; base; base= base->next) {
 346                 if( (base->lay & layer) && base->object->pd) {
 347                         if(base->object->pd->forcefield) retval |= USES_FIELD;
 348                         if(base->object->pd->deflect) retval |= USES_DEFLECT;
 349                 }
 350         }
 351         return retval;
 352 }
 353
 354 static void softbody_calc_forces(Object *ob, float forcetime)
 355 {
 356 /* rule we never alter free variables :bp->vec bp->pos in here !
 357  * this will ruin adaptive stepsize AKA heun! (BM)
 358  */
 359         SoftBody *sb= ob->soft; // is supposed to be there
 360         BodyPoint  *bp;
 361         BodyPoint *bproot;
 362         BodySpring *bs;
 363         float iks, ks, kd, gravity, actspringlen, forcefactor, sd[3],
 364         fieldfactor = 1000.0f,
 365         windfactor  = 250.0f;
 366         int a, b, do_effector;
 367
 368         /* clear forces */
 369         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 370                 bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
 371         }
 372
 373         gravity = sb->grav * sb_grav_force_scale(ob);
 374         /* check! */
 375         do_effector= is_there_deflection(ob->lay);
 376         iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
 377         bproot= sb->bpoint; /* need this for proper spring addressing */
 378
 379         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 380                 if(bp->goal < SOFTGOALSNAP){ // ommit this bp when i snaps
 381                         float auxvect[3]; // aux unit vector
 382                         float velgoal[3];
 383                         float absvel =0, projvel= 0;
 384
 385                         /* do goal stuff */
 386                         if(ob->softflag & OB_SB_GOAL) {
 387                                 /* true elastic goal */
 388                                 VecSubf(auxvect,bp->origT,bp->pos);
 389                                 ks  = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
 390                                 bp->force[0]= ks*(auxvect[0]);
 391                                 bp->force[1]= ks*(auxvect[1]);
 392                                 bp->force[2]= ks*(auxvect[2]);
 393                                 /* calulate damping forces generated by goals*/
 394                                 VecSubf(velgoal,bp->origS, bp->origE);
 395                                 kd =  sb->goalfrict * sb_fric_force_scale(ob) ;
 396
 397                                 if (forcetime > 0.0 ) { // make sure friction does not become rocket motor on time reversal
 398                                         bp->force[0]-= kd * (velgoal[0] + bp->vec[0]);
 399                                         bp->force[1]-= kd * (velgoal[1] + bp->vec[1]);
 400                                         bp->force[2]-= kd * (velgoal[2] + bp->vec[2]);
 401                                 }
 402                                 else {
 403                                         bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
 404                                         bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
 405                                         bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
 406                                 }
 407                         }
 408                         /* done goal stuff */
 409
 410
 411                         /* gravitation */
 412                         bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
 413
 414                         /* particle field & vortex */
 415                         if(do_effector & USES_FIELD) {
 416                                 float force[3]= {0.0f, 0.0f, 0.0f};
 417                                 float speed[3]= {0.0f, 0.0f, 0.0f};
 418                                 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); // just for calling functio once
 419
 420                                 pdDoEffector(bp->pos, force, speed, (float)G.scene->r.cfra, ob->lay,PE_WIND_AS_SPEED);
 421                                 // note: now we have wind as motion of media, so we can do anisotropic stuff here,
 422                                 // if we had vertex normals here(BM)
 423                                 /* apply forcefield*/
 424                                 VecMulf(force,fieldfactor* eval_sb_fric_force_scale);
 425                                 VECADD(bp->force, bp->force, force);
 426
 427                                 /* friction in moving media */
 428                                 kd= sb->mediafrict* eval_sb_fric_force_scale;
 429                                 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
 430                                 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
 431                                 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
 432                                 /* now we'll have nice centrifugal effect for vortex */
 433
 434                         }
 435                         else {
 436                                 /* friction in media (not) moving*/
 437                                 kd= sb->mediafrict* sb_fric_force_scale(ob);
 438                                 /* assume it to be proportional to actual velocity */
 439                                 bp->force[0]-= bp->vec[0]*kd;
 440                                 bp->force[1]-= bp->vec[1]*kd;
 441                                 bp->force[2]-= bp->vec[2]*kd;
 442                                 /* friction in media done */
 443                         }
 444
 445                         /*other forces*/
 446                         /* this is the place where other forces can be added
 447                         yes, constraints and collision stuff should go here too (read baraff papers on that!)
 448                         */
 449                         /* try to match moving collision targets */
 450                         /* master switch to turn collision off (BM)*/
 451                         //if(0) {
 452                         if(do_effector & USES_DEFLECT) {
 453                                 /*sorry for decl. here i'll move 'em up when WIP is done (BM) */
 454                                 float defforce[3];
 455                                 float collisionpos[3],facenormal[3];
 456                                 float cf = 1.0f;
 457                                 float bounce = 0.5f;
 458                                 kd = 1.0f;
 459                                 defforce[0] = 0.0f;
 460                                 defforce[1] = 0.0f;
 461                                 defforce[2] = 0.0f;
 462
 463                                 if (sb_deflect_face(ob,bp->pos, bp->pos, collisionpos, facenormal,defforce,&cf,&bounce)){
 464                                         bp->force[0] += defforce[0]*kd;
 465                                         bp->force[1] += defforce[1]*kd;
 466                                         bp->force[2] += defforce[2]*kd;
 467                                         bp->contactfrict = cf;
 468                                 }
 469                                 else{
 470                                         bp->contactfrict = 0.0f;
 471                                 }
 472
 473                         }
 474
 475                         /*other forces done*/
 476                         /* nice things could be done with anisotropic friction
 477                         like wind/air resistance in normal direction
 478                         --> having a piece of cloth sailing down
 479                         but this needs to have a *valid* vertex normal
 480                         *valid* means to be calulated on time axis
 481                         hrms .. may be a rough one could be used as well .. let's see
 482                         */
 483
 484                         if(ob->softflag & OB_SB_EDGES) {
 485                                 if (sb->bspring){ // spring list exists at all ?
 486                                         for(b=bp->nofsprings;b>0;b--){
 487                                                 bs = sb->bspring + bp->springs[b-1];
 488                                                 if (( (sb->totpoint-a) == bs->v1) ){
 489                                                         actspringlen= VecLenf( (bproot+bs->v2)->pos, bp->pos);
 490                                                         VecSubf(sd,(bproot+bs->v2)->pos, bp->pos);
 491                                                         Normalise(sd);
 492
 493                                                         // friction stuff V1
 494                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v2)->vec);
 495                                                         kd = sb->infrict * sb_fric_force_scale(ob);
 496                                                         absvel  = Normalise(velgoal);
 497                                                         projvel = ABS(Inpf(sd,velgoal));
 498                                                         kd *= absvel * projvel;
 499                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
 500
 501                                                         if(bs->len > 0.0) /* check for degenerated springs */
 502                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
 503                                                         else
 504                                                                 forcefactor = actspringlen * iks;
 505
 506                                                         Vec3PlusStVec(bp->force,-forcefactor,sd);
 507
 508                                                 }
 509
 510                                                 if (( (sb->totpoint-a) == bs->v2) ){
 511                                                         actspringlen= VecLenf( (bproot+bs->v1)->pos, bp->pos);
 512                                                         VecSubf(sd,bp->pos,(bproot+bs->v1)->pos);
 513                                                         Normalise(sd);
 514
 515                                                         // friction stuff V2
 516                                                         VecSubf(velgoal,bp->vec,(bproot+bs->v1)->vec);
 517                                                         kd = sb->infrict * sb_fric_force_scale(ob);
 518                                                         absvel  = Normalise(velgoal);
 519                                                         projvel = ABS(Inpf(sd,velgoal));
 520                                                         kd *= absvel * projvel;
 521                                                         Vec3PlusStVec(bp->force,-kd,velgoal);
 522
 523                                                         if(bs->len > 0.0)
 524                                                                 forcefactor = (bs->len - actspringlen)/bs->len * iks;
 525                                                         else
 526                                                                 forcefactor = actspringlen * iks;
 527                                                         Vec3PlusStVec(bp->force,+forcefactor,sd);
 528                                                 }
 529                                         }/* loop springs */
 530                                 }/* existing spring list */
 531                         }/*any edges*/
 532                 }/*omit on snap */
 533         }/*loop all bp's*/
 534 }
 535
 536 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err)
 537 {
 538         /* time evolution */
 539         /* actually does an explicit euler step mode == 0 */
 540         /* or heun ~ 2nd order runge-kutta steps, mode 1,2 */
 541         SoftBody *sb= ob->soft; // is supposed to be there
 542         BodyPoint *bp;
 543         float dx[3],dv[3];
 544         float timeovermass;
 545         float maxerr = 0.0;
 546         int a, do_effector;
 547
 548     forcetime *= sb_time_scale(ob);
 549         /* check! */
 550         do_effector= is_there_deflection(ob->lay);
 551
 552         // claim a minimum mass for vertex
 553         if (sb->nodemass > 0.09999f) timeovermass = forcetime/sb->nodemass;
 554         else timeovermass = forcetime/0.09999f;
 555
 556         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 557                 if(bp->goal < SOFTGOALSNAP){
 558
 559                         /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces  + more forces*/
 560                         /* the ( ... )' operator denotes derivate respective time */
 561                         /* the euler step for velocity then becomes */
 562                         /* v(t + dt) = v(t) + a(t) * dt */
 563                         bp->force[0]*= timeovermass; /* individual mass of node here */
 564                         bp->force[1]*= timeovermass;
 565                         bp->force[2]*= timeovermass;
 566                         /* some nasty if's to have heun in here too */
 567                         VECCOPY(dv,bp->force);
 568                         if (mode == 1){
 569                                 VECCOPY(bp->prevvec, bp->vec);
 570                                 VECCOPY(bp->prevdv, dv);
 571                         }
 572                         if (mode ==2){
 573                                 /* be optimistic and execute step */
 574                                 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]);
 575                                 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]);
 576                                 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]);
 577                                 /* compare euler to heun to estimate error for step sizing */
 578                                 maxerr = MAX2(maxerr,ABS(dv[0] - bp->prevdv[0]));
 579                                 maxerr = MAX2(maxerr,ABS(dv[1] - bp->prevdv[1]));
 580                                 maxerr = MAX2(maxerr,ABS(dv[2] - bp->prevdv[2]));
 581                         }
 582                         else {VECADD(bp->vec, bp->vec, bp->force);}
 583
 584                         /* so here is (x)'= v(elocity) */
 585                         /* the euler step for location then becomes */
 586                         /* x(t + dt) = x(t) + v(t) * dt */
 587
 588                         VECCOPY(dx,bp->vec);
 589                         dx[0]*=forcetime ;
 590                         dx[1]*=forcetime ;
 591                         dx[2]*=forcetime ;
 592
 593                         /* again some nasty if's to have heun in here too */
 594                         if (mode ==1){
 595                                 VECCOPY(bp->prevpos,bp->pos);
 596                                 VECCOPY(bp->prevdx ,dx);
 597                         }
 598
 599                         if (mode ==2){
 600                                 bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]);
 601                                 bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]);
 602                                 bp->pos[2] = bp->prevpos[2] + 0.5f* ( dx[2] + bp->prevdx[2]);
 603                                 maxerr = MAX2(maxerr,ABS(dx[0] - bp->prevdx[0]));
 604                                 maxerr = MAX2(maxerr,ABS(dx[1] - bp->prevdx[1]));
 605                                 maxerr = MAX2(maxerr,ABS(dx[2] - bp->prevdx[2]));
 606 /* kind of hack .. while inside collision target .. make movement more *viscous* */
 607                                 if (bp->contactfrict > 0.0f){
 608                                         bp->vec[0] *= (1.0 - bp->contactfrict);
 609                                         bp->vec[1] *= (1.0 - bp->contactfrict);
 610                                         bp->vec[2] *= (1.0 - bp->contactfrict);
 611                                 }
 612                         }
 613                         else { VECADD(bp->pos, bp->pos, dx);}
 614 // experimental particle collision suff was here .. just to help HOS on next merge (BM)
 615                 }//snap
 616         } //for
 617         if (err){ /* so step size will be controlled by biggest difference in slope */
 618                 *err = maxerr;
 619         }
 620 }
 621
 622 /* used by heun when it overshoots */
 623 static void softbody_restore_prev_step(Object *ob)
 624 {
 625         SoftBody *sb= ob->soft; // is supposed to be there
 626         BodyPoint *bp;
 627         int a;
 628
 629         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 630                 VECCOPY(bp->vec, bp->prevvec);
 631                 VECCOPY(bp->pos, bp->prevpos);
 632         }
 633 }
 634
 635
 636 static void softbody_apply_goalsnap(Object *ob)
 637 {
 638         SoftBody *sb= ob->soft; // is supposed to be there
 639         BodyPoint *bp;
 640         int a;
 641
 642         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 643                 if (bp->goal >= SOFTGOALSNAP){
 644                         VECCOPY(bp->prevpos,bp->pos);
 645                         VECCOPY(bp->pos,bp->origT);
 646                 }
 647         }
 648 }
 649
 650 /* unused */
 651 #if 0
 652 static void softbody_force_goal(Object *ob)
 653 {
 654         SoftBody *sb= ob->soft; // is supposed to be there
 655         BodyPoint *bp;
 656         int a;
 657
 658         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 659                 VECCOPY(bp->pos,bp->origT);
 660                 bp->vec[0] = bp->origE[0] - bp->origS[0];
 661                 bp->vec[1] = bp->origE[1] - bp->origS[1];
 662                 bp->vec[2] = bp->origE[2] - bp->origS[2];
 663         }
 664 }
 665 #endif
 666
 667 /* expects full initialized softbody */
 668 static void interpolate_exciter(Object *ob, int timescale, int time)
 669 {
 670         SoftBody *sb= ob->soft;
 671         BodyPoint *bp;
 672         float f;
 673         int a;
 674
 675         // note: i removed Mesh usage here, softbody should remain generic! (ton)
 676
 677         f = (float)time/(float)timescale;
 678
 679         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
 680                 bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]);
 681                 bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]);
 682                 bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]);
 683                 if (bp->goal >= SOFTGOALSNAP){
 684                         bp->vec[0] = bp->origE[0] - bp->origS[0];
 685                         bp->vec[1] = bp->origE[1] - bp->origS[1];
 686                         bp->vec[2] = bp->origE[2] - bp->origS[2];
 687                 }
 688         }
 689
 690         if(ob->softflag & OB_SB_EDGES) {
 691                 /* hrms .. do springs alter their lenght ?
 692                 bs= ob->soft->bspring;
 693                 bp= ob->soft->bpoint;
 694                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, bs++) {
 695                         bs->len= VecLenf( (bp+bs->v1)->origT, (bp+bs->v2)->origT);
 696                 }
 697                 */
 698         }
 699 }
 700
 701
 702 /* ************ convertors ********** */
 703
 704 /*  for each object type we need;
 705     - xxxx_to_softbody(Object *ob)      : a full (new) copy
 706         - xxxx_update_softbody(Object *ob)  : update refreshes current positions
 707     - softbody_to_xxxx(Object *ob)      : after simulation, copy vertex locations back
 708 */
 709
 710 /* helper  call */
 711 static int object_has_edges(Object *ob)
 712 {
 713         if(ob->type==OB_MESH) {
 714                 Mesh *me= ob->data;
 715                 if(me->medge) return 1;
 716         }
 717         else if(ob->type==OB_LATTICE) {
 718                 ;
 719         }
 720
 721         return 0;
 722 }
 723
 724 /* helper  call */
 725 static void set_body_point(Object *ob, BodyPoint *bp, float *vec)
 726 {
 727
 728         VECCOPY(bp->pos, vec);
 729         Mat4MulVecfl(ob->obmat, bp->pos);  // yep, sofbody is global coords
 730         VECCOPY(bp->origS, bp->pos);
 731         VECCOPY(bp->origE, bp->pos);
 732         VECCOPY(bp->origT, bp->pos);
 733
 734         bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0;
 735         bp->weight= 1.0;
 736         bp->goal= ob->soft->defgoal;
 737
 738         bp->nofsprings= 0;
 739         bp->springs= NULL;
 740         bp->contactfrict = 0.0f;
 741 }
 742
 743
 744 /* copy original (new) situation in softbody, as result of matrices or deform */
 745 /* is assumed to enter function with ob->soft, but can be without points */
 746 static void mesh_update_softbody(Object *ob)
 747 {
 748         Mesh *me= ob->data;
 749         MVert *mvert= me->mvert;
 750 /*      MEdge *medge= me->medge;  */ /*unused*/
 751         BodyPoint *bp;
 752         int a;
 753
 754         /* possible after a file read... */
 755         if(ob->soft->totpoint!=me->totvert) sbObjectToSoftbody(ob);
 756
 757         if(me->totvert) {
 758
 759                 bp= ob->soft->bpoint;
 760                 for(a=0; a<me->totvert; a++, mvert++, bp++) {
 761                         VECCOPY(bp->origS, bp->origE);
 762                         VECCOPY(bp->origE, mvert->co);
 763                         Mat4MulVecfl(ob->obmat, bp->origE);
 764                         VECCOPY(bp->origT, bp->origE);
 765                 }
 766
 767                 if(ob->softflag & OB_SB_EDGES) {
 768
 769                         /* happens when in UI edges was set */
 770                         if(ob->soft->bspring==NULL)
 771                                 if(object_has_edges(ob)) sbObjectToSoftbody(ob);
 772
 773                         /* hrms .. do springs alter their lenght ? (yes, mesh keys would (ton))
 774                         if(medge) {
 775                                 bs= ob->soft->bspring;
 776                                 bp= ob->soft->bpoint;
 777                                 for(a=0; (a<me->totedge && a < ob->soft->totspring ); a++, medge++, bs++) {
 778                                         bs->len= VecLenf( (bp+bs->v1)->origE, (bp+bs->v2)->origE);
 779                                 }
 780                         }
 781                         */
 782                 }
 783         }
 784 }
 785
 786
 787 static void get_scalar_from_vertexgroup(Object *ob, int vertID, short groupindex, float *target)
 788 /* result 0 on success, else indicates error number
 789 -- kind of *inverse* result defintion,
 790 -- but this way we can signal error condition to caller
 791 -- and yes this function must not be here but in a *vertex group module*
 792 */
 793 {
 794         MDeformVert *dv;
 795         int i;
 796
 797         /* spot the vert in deform vert list at mesh */
 798         if(ob->type==OB_MESH) {
 799                 if (((Mesh *)ob->data)->dvert) {
 800                         dv = ((Mesh*)ob->data)->dvert + vertID;
 801                         /* Lets see if this vert is in the weight group */
 802                         for (i=0; i<dv->totweight; i++){
 803                                 if (dv->dw[i].def_nr == groupindex){
 804                                         *target= dv->dw[i].weight; /* got it ! */
 805                                         break;
 806                                 }
 807                         }
 808                 }
 809         }
 810 }
 811
 812 /* makes totally fresh start situation */
 813 static void mesh_to_softbody(Object *ob)
 814 {
 815         SoftBody *sb;
 816         Mesh *me= ob->data;
 817         MVert *mvert= me->mvert;
 818         MEdge *medge= me->medge;
 819         BodyPoint *bp;
 820         BodySpring *bs;
 821         float goalfac;
 822         int a, totedge;
 823
 824         if (ob->softflag & OB_SB_EDGES) totedge= me->totedge;
 825         else totedge= 0;
 826
 827         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
 828         renew_softbody(ob, me->totvert, totedge);
 829
 830         /* we always make body points */
 831         sb= ob->soft;
 832         bp= sb->bpoint;
 833         goalfac= ABS(sb->maxgoal - sb->mingoal);
 834
 835         for(a=me->totvert; a>0; a--, mvert++, bp++) {
 836
 837                 set_body_point(ob, bp, mvert->co);
 838
 839                 /* get scalar values needed  *per vertex* from vertex group functions,
 840                 so we can *paint* them nicly ..
 841                 they are normalized [0.0..1.0] so may be we need amplitude for scale
 842                 which can be done by caller but still .. i'd like it to go this way
 843                 */
 844
 845                 if(sb->vertgroup) {
 846                         get_scalar_from_vertexgroup(ob, me->totvert - a, sb->vertgroup-1, &bp->goal);
 847                         // do this always, regardless successfull read from vertex group
 848                         bp->goal= sb->mingoal + bp->goal*goalfac;
 849                 }
 850                 /* a little ad hoc changing the goal control to be less *sharp* */
 851                 bp->goal = (float)pow(bp->goal, 4.0f);
 852
 853                 /* to proove the concept
 854                 this would enable per vertex *mass painting*
 855                 strcpy(name,"SOFTMASS");
 856                 error = get_scalar_from_named_vertexgroup(ob,name,me->totvert - a,&temp);
 857                 if (!error) bp->mass = temp * ob->rangeofmass;
 858                 */
 859         }
 860
 861         /* but we only optionally add body edge springs */
 862         if (ob->softflag & OB_SB_EDGES) {
 863                 if(medge) {
 864                         bs= sb->bspring;
 865                         bp= sb->bpoint;
 866                         for(a=me->totedge; a>0; a--, medge++, bs++) {
 867                                 bs->v1= medge->v1;
 868                                 bs->v2= medge->v2;
 869                                 bs->strength= 1.0;
 870                                 bs->len= VecLenf( (bp+bs->v1)->origS, (bp+bs->v2)->origS);
 871                         }
 872
 873
 874                         /* insert *diagonal* springs in quads if desired */
 875                         if (ob->softflag & OB_SB_QUADS) {
 876                                 add_mesh_quad_diag_springs(ob);
 877                         }
 878
 879                         build_bps_springlist(ob); /* big mesh optimization */
 880                 }
 881         }
 882
 883 }
 884
 885 /* copies current sofbody position in mesh, so do this within modifier stacks! */
 886 static void softbody_to_mesh(Object *ob)
 887 {
 888         Mesh *me= ob->data;
 889         MVert *mvert;
 890         BodyPoint *bp;
 891         int a;
 892
 893         bp= ob->soft->bpoint;
 894         mvert= me->mvert;
 895         for(a=me->totvert; a>0; a--, mvert++, bp++) {
 896                 VECCOPY(mvert->co, bp->pos);
 897                 Mat4MulVecfl(ob->imat, mvert->co);      // softbody is in global coords
 898         }
 899
 900 }
 901
 902 /* makes totally fresh start situation */
 903 static void lattice_to_softbody(Object *ob)
 904 {
 905         SoftBody *sb;
 906         Lattice *lt= ob->data;
 907         BodyPoint *bop;
 908         BPoint *bp;
 909         int a, totvert;
 910
 911         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
 912
 913         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
 914         renew_softbody(ob, totvert, 0);
 915
 916         /* we always make body points */
 917         sb= ob->soft;
 918
 919         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
 920                 set_body_point(ob, bop, bp->vec);
 921         }
 922 }
 923
 924 /* copies current sofbody position */
 925 static void softbody_to_lattice(Object *ob)
 926 {
 927         SoftBody *sb;
 928         Lattice *lt= ob->data;
 929         BodyPoint *bop;
 930         BPoint *bp;
 931         int a, totvert;
 932
 933         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
 934         sb= ob->soft;
 935
 936         for(a= totvert, bp= lt->def, bop= sb->bpoint; a>0; a--, bp++, bop++) {
 937                 VECCOPY(bp->vec, bop->pos);
 938                 Mat4MulVecfl(ob->imat, bp->vec);        // softbody is in global coords
 939         }
 940 }
 941
 942 /* copy original (new) situation in softbody, as result of matrices or deform */
 943 /* is assumed to enter function with ob->soft, but can be without points */
 944 static void lattice_update_softbody(Object *ob)
 945 {
 946         Lattice *lt= ob->data;
 947         BodyPoint *bop;
 948         BPoint *bp;
 949         int a, totvert;
 950
 951         totvert= lt->pntsu*lt->pntsv*lt->pntsw;
 952
 953         /* possible after a file read... */
 954         if(ob->soft->totpoint!=totvert) sbObjectToSoftbody(ob);
 955
 956         for(a= totvert, bp= lt->def, bop= ob->soft->bpoint; a>0; a--, bp++, bop++) {
 957                 VECCOPY(bop->origS, bop->origE);
 958                 VECCOPY(bop->origE, bp->vec);
 959                 Mat4MulVecfl(ob->obmat, bop->origE);
 960                 VECCOPY(bop->origT, bop->origE);
 961         }
 962 }
 963
 964
 965 /* copies softbody result back in object */
 966 /* only used in sbObjectStep() */
 967 static void softbody_to_object(Object *ob)
 968 {
 969
 970         if(ob->soft==NULL) return;
 971
 972         /* inverse matrix is not uptodate... */
 973         Mat4Invert(ob->imat, ob->obmat);
 974
 975         switch(ob->type) {
 976         case OB_MESH:
 977                 softbody_to_mesh(ob);
 978                 break;
 979         case OB_LATTICE:
 980                 softbody_to_lattice(ob);
 981                 break;
 982         }
 983 }
 984
 985 /* copy original (new) situation in softbody, as result of matrices or deform */
 986 /* used in sbObjectStep() and sbObjectReset() */
 987 /* assumes to have ob->soft, but can be entered without points */
 988 static void object_update_softbody(Object *ob)
 989 {
 990
 991         switch(ob->type) {
 992         case OB_MESH:
 993                 mesh_update_softbody(ob);
 994                 break;
 995         case OB_LATTICE:
 996                 lattice_update_softbody(ob);
 997                 break;
 998         }
 999
1000 }
1001
1002
1003
1004 /* ************ Object level, exported functions *************** */
1005
1006 /* allocates and initializes general main data */
1007 SoftBody *sbNew(void)
1008 {
1009         SoftBody *sb;
1010
1011         sb= MEM_callocN(sizeof(SoftBody), "softbody");
1012
1013         sb->mediafrict= 0.5;
1014         sb->nodemass= 1.0;
1015         sb->grav= 0.0;
1016         sb->physics_speed= 1.0;
1017         sb->rklimit= 0.1;
1018
1019         sb->goalspring= 0.5;
1020         sb->goalfrict= 0.0;
1021         sb->mingoal= 0.0;
1022         sb->maxgoal= 1.0;
1023         sb->defgoal= 0.7;
1024
1025         sb->inspring= 0.5;
1026         sb->infrict= 0.5;
1027
1028         return sb;
1029 }
1030
1031 /* frees all */
1032 void sbFree(SoftBody *sb)
1033 {
1034         free_softbody_intern(sb);
1035         MEM_freeN(sb);
1036 }
1037
1038
1039 /* makes totally fresh start situation */
1040 void sbObjectToSoftbody(Object *ob)
1041 {
1042
1043         switch(ob->type) {
1044         case OB_MESH:
1045                 mesh_to_softbody(ob);
1046                 break;
1047         case OB_LATTICE:
1048                 lattice_to_softbody(ob);
1049                 break;
1050         }
1051
1052         if(ob->soft) ob->soft->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1053         ob->softflag &= ~OB_SB_REDO;
1054 }
1055
1056 /* reset all motion */
1057 void sbObjectReset(Object *ob)
1058 {
1059         SoftBody *sb= ob->soft;
1060         BodyPoint *bp;
1061         int a;
1062
1063         if(sb==NULL) return;
1064         sb->ctime= bsystem_time(ob, NULL, (float)G.scene->r.cfra, 0.0);
1065
1066         object_update_softbody(ob);
1067
1068         for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
1069                 // origS is previous timestep
1070                 VECCOPY(bp->origS, bp->origE);
1071                 VECCOPY(bp->pos, bp->origE);
1072                 VECCOPY(bp->origT, bp->origE);
1073                 bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f;
1074
1075                 // no idea about the Heun stuff! (ton)
1076                 VECCOPY(bp->prevpos, bp->pos);
1077                 VECCOPY(bp->prevvec, bp->vec);
1078                 VECCOPY(bp->prevdx, bp->vec);
1079                 VECCOPY(bp->prevdv, bp->vec);
1080         }
1081 }
1082
1083
1084 /* simulates one step. framenr is in frames */
1085 /* copies result back to object, displist */
1086 void sbObjectStep(Object *ob, float framenr)
1087 {
1088         SoftBody *sb;
1089         Base *base;
1090         float dtime;
1091         int timescale,t;
1092         float ctime, forcetime;
1093         float err;
1094
1095         /* this variable is set while transform(). with lattices also having a softbody,
1096            it calls lattice_modifier() all the time... has no displist yet. Is temporal
1097            hack which should be resolved with proper depgraph usage + storage of deformed
1098            vertices in lattice (ton) */
1099         if(G.moving) return;
1100
1101         /* remake softbody if: */
1102         if( (ob->softflag & OB_SB_REDO) ||              // signal after weightpainting
1103                 (ob->soft==NULL) ||                                     // just to be nice we allow full init
1104                 (ob->soft->totpoint==0) )                       // after reading new file, or acceptable as signal to refresh
1105                         sbObjectToSoftbody(ob);
1106
1107         sb= ob->soft;
1108
1109         /* still no points? go away */
1110         if(sb->totpoint==0) return;
1111
1112         /* reset deflector cache, sumohandle is free, but its still sorta abuse... (ton) */
1113         for(base= G.scene->base.first; base; base= base->next) {
1114                 base->object->sumohandle= NULL;
1115         }
1116
1117         /* checking time: */
1118         ctime= bsystem_time(ob, NULL, framenr, 0.0);
1119         dtime= ctime - sb->ctime;
1120                 // bail out for negative or for large steps
1121         if(dtime<0.0 || dtime >= 9.9*G.scene->r.framelen) { // note: what is G.scene->r.framelen for ? (BM)
1122                 sbObjectReset(ob);
1123                 return;
1124         }
1125
1126         /* the simulator */
1127
1128         if(dtime > 0.0) {       // note: what does this mean now? (ton)
1129                 //answer (BM) :
1130                 //dtime is still in [frames]
1131                 //we made sure dtime is >= 0.0
1132                 //but still need to handle dtime == 0.0 -> just return sb as is, just to be nice
1133                 object_update_softbody(ob);
1134
1135                 if (TRUE) {     // RSOL1 always true now (ton)
1136                         /* special case of 2nd order Runge-Kutta type AKA Heun */
1137                         float timedone =0.0; // how far did we get without violating error condition
1138                         /* loops = counter for emergency brake
1139                          * we don't want to lock up the system if physics fail
1140                          */
1141                         int loops =0 ;
1142                         SoftHeunTol = sb->rklimit; // humm .. this should be calculated from sb parameters and sizes
1143
1144                         forcetime = dtime; /* hope for integrating in one step */
1145                         while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) )
1146                         {
1147                                 if (ABS(dtime) > 9.0 ){
1148                                         if(G.f & G_DEBUG) printf("SB_STEPSIZE \n");
1149                                         break; // sorry but i must assume goal movement can't be interpolated any more
1150                                 }
1151                                 //set goals in time
1152                                 interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime)));
1153                                 // do predictive euler step
1154                                 softbody_calc_forces(ob, forcetime);
1155                                 softbody_apply_forces(ob, forcetime, 1, NULL);
1156                                 // crop new slope values to do averaged slope step
1157                                 softbody_calc_forces(ob, forcetime);
1158                                 softbody_apply_forces(ob, forcetime, 2, &err);
1159                                 softbody_apply_goalsnap(ob);
1160
1161                                 if (err > SoftHeunTol){ // error needs to be scaled to some quantity
1162                                         softbody_restore_prev_step(ob);
1163                                         forcetime /= 2.0;
1164                                 }
1165                                 else {
1166                                         float newtime = forcetime * 1.1f; // hope for 1.1 times better conditions in next step
1167
1168                                         if (err > SoftHeunTol/2.0){ // stay with this stepsize unless err really small
1169                                                 newtime = forcetime;
1170                                         }
1171                                         timedone += forcetime;
1172                                         if (forcetime > 0.0)
1173                                                 forcetime = MIN2(dtime - timedone,newtime);
1174                                         else
1175                                                 forcetime = MAX2(dtime - timedone,newtime);
1176                                 }
1177                                 loops++;
1178                         }
1179                         // move snapped to final position
1180                         interpolate_exciter(ob, 2, 2);
1181                         softbody_apply_goalsnap(ob);
1182
1183                         if(G.f & G_DEBUG) {
1184                                 if (loops > HEUNWARNLIMIT) /* monitor high loop counts say 1000 after testing */
1185                                         printf("%d heun integration loops/frame \n",loops);
1186                         }
1187                 }
1188                 else{
1189                         /* do brute force explicit euler */
1190                         /* inner intagration loop */
1191                         /* */
1192                         // loop n times so that n*h = duration of one frame := 1
1193                         // x(t+h) = x(t) + h*v(t);
1194                         // v(t+h) = v(t) + h*f(x(t),t);
1195                         timescale = (int)(sb->rklimit * ABS(dtime));
1196                         for(t=1 ; t <= timescale; t++) {
1197                                 if (ABS(dtime) > 15 ) break;
1198
1199                                 /* the *goal* mesh must use the n*h timing too !
1200                                 use *cheap* linear intepolation for that  */
1201                                 interpolate_exciter(ob,timescale,t);
1202                                 if (timescale > 0 ) {
1203                                         forcetime = dtime/timescale;
1204
1205                                         /* does not fit the concept sloving ODEs :) */
1206                                         /*                      softbody_apply_goal(ob,forcetime );  */
1207
1208                                         /* explicit Euler integration */
1209                                         /* we are not controling a nuclear power plant!
1210                                         so rought *almost* physical behaviour is acceptable.
1211                                         in cases of *mild* stiffnes cranking up timscale -> decreasing stepsize *h*
1212                                         avoids instability      */
1213                                         softbody_calc_forces(ob,forcetime);
1214                                         softbody_apply_forces(ob,forcetime,0, NULL);
1215                                         softbody_apply_goalsnap(ob);
1216
1217                                         //      if (0){
1218                                         /* ok here comes the überhammer
1219                                         use a semi implicit euler integration to tackle *all* stiff conditions
1220                                         but i doubt the cost/benifit holds for most of the cases
1221                                         -- to be coded*/
1222                                         //      }
1223
1224                                 }
1225                         }
1226                 }
1227
1228                 /* and apply to vertices */
1229                  softbody_to_object(ob);
1230
1231                 sb->ctime= ctime;
1232         } // if(ABS(dtime) > 0.0)
1233         else {
1234                 // rule : you have asked for the current state of the softobject
1235                 // since dtime= ctime - ob->soft->ctime== 0.0;
1236                 // and we were not notifified about any other time changes
1237                 // so here it is !
1238                 softbody_to_object(ob);
1239         }
1240
1241         /* reset deflector cache */
1242         for(base= G.scene->base.first; base; base= base->next) {
1243                 if(base->object->sumohandle) {
1244                         MEM_freeN(base->object->sumohandle);
1245                         base->object->sumohandle= NULL;
1246                 }
1247         }
1248
1249 }
1250