Cleanup: Remove more #if 0 blocks
[blender.git] / source / blender / blenkernel / intern / softbody.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) Blender Foundation
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): none yet.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenkernel/intern/softbody.c
29  *  \ingroup bke
30  */
31
32
33 /**
34  * variables on the UI for now
35  * <pre>
36  * float mediafrict;  friction to env
37  * float nodemass;        softbody mass of *vertex*
38  * float grav;        softbody amount of gravitaion to apply
39  *
40  * float goalspring;  softbody goal springs
41  * float goalfrict;   softbody goal springs friction
42  * float mingoal;     quick limits for goal
43  * float maxgoal;
44  *
45  * float inspring;        softbody inner springs
46  * float infrict;     softbody inner springs friction
47  * </pre>
48  */
49
50
51 #include <math.h>
52 #include <stdlib.h>
53 #include <string.h>
54
55 #include "MEM_guardedalloc.h"
56
57 /* types */
58 #include "DNA_collection_types.h"
59 #include "DNA_curve_types.h"
60 #include "DNA_lattice_types.h"
61 #include "DNA_meshdata_types.h"
62 #include "DNA_mesh_types.h"
63 #include "DNA_object_types.h"
64 #include "DNA_scene_types.h"
65
66 #include "BLI_math.h"
67 #include "BLI_utildefines.h"
68 #include "BLI_listbase.h"
69 #include "BLI_ghash.h"
70 #include "BLI_threads.h"
71
72 #include "BKE_collection.h"
73 #include "BKE_collision.h"
74 #include "BKE_curve.h"
75 #include "BKE_effect.h"
76 #include "BKE_global.h"
77 #include "BKE_layer.h"
78 #include "BKE_modifier.h"
79 #include "BKE_softbody.h"
80 #include "BKE_pointcache.h"
81 #include "BKE_deform.h"
82 #include "BKE_mesh.h"
83 #include "BKE_scene.h"
84
85 #include "DEG_depsgraph.h"
86 #include "DEG_depsgraph_query.h"
87
88 #include  "PIL_time.h"
89
90 /* callbacks for errors and interrupts and some goo */
91 static int (*SB_localInterruptCallBack)(void) = NULL;
92
93
94 /* ********** soft body engine ******* */
95
96 typedef enum {SB_EDGE=1, SB_BEND=2, SB_STIFFQUAD=3, SB_HANDLE=4} type_spring;
97
98 typedef struct BodySpring {
99         int v1, v2;
100         float len, cf, load;
101         float ext_force[3]; /* edges colliding and sailing */
102         type_spring springtype;
103         short flag;
104 } BodySpring;
105
106 typedef struct BodyFace {
107         int v1, v2, v3;
108         float ext_force[3]; /* faces colliding */
109         short flag;
110 } BodyFace;
111
112 typedef struct ReferenceVert {
113         float pos[3]; /* position relative to com */
114         float mass;   /* node mass */
115 } ReferenceVert;
116
117 typedef struct ReferenceState {
118         float com[3]; /* center of mass*/
119         ReferenceVert *ivert; /* list of initial values */
120 } ReferenceState;
121
122
123 /*private scratch pad for caching and other data only needed when alive*/
124 typedef struct SBScratch {
125         GHash *colliderhash;
126         short needstobuildcollider;
127         short flag;
128         BodyFace *bodyface;
129         int totface;
130         float aabbmin[3], aabbmax[3];
131         ReferenceState Ref;
132 } SBScratch;
133
134 typedef struct  SB_thread_context {
135                 Scene *scene;
136                 Object *ob;
137                 float forcetime;
138                 float timenow;
139                 int ifirst;
140                 int ilast;
141                 ListBase *effectors;
142                 int do_deflector;
143                 float fieldfactor;
144                 float windfactor;
145                 int nr;
146                 int tot;
147 } SB_thread_context;
148
149 #define MID_PRESERVE 1
150
151 #define SOFTGOALSNAP  0.999f
152 /* if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
153  * removes *unnecessary* stiffness from ODE system
154  */
155 #define HEUNWARNLIMIT 1 /* 500 would be fine i think for detecting severe *stiff* stuff */
156
157
158 #define BSF_INTERSECT   1 /* edge intersects collider face */
159
160 /* private definitions for bodypoint states */
161 #define SBF_DOFUZZY          1 /* Bodypoint do fuzzy */
162 #define SBF_OUTOFCOLLISION   2 /* Bodypoint does not collide  */
163
164
165 #define BFF_INTERSECT   1 /* collider edge   intrudes face */
166 #define BFF_CLOSEVERT   2 /* collider vertex repulses face */
167
168
169 static float SoftHeunTol = 1.0f; /* humm .. this should be calculated from sb parameters and sizes */
170
171 /* local prototypes */
172 static void free_softbody_intern(SoftBody *sb);
173
174 /*+++ frame based timing +++*/
175
176 /*physical unit of force is [kg * m / sec^2]*/
177
178 static float sb_grav_force_scale(Object *UNUSED(ob))
179 /* since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm
180  * put it to a function here, so we can add user options later without touching simulation code
181  */
182 {
183         return (0.001f);
184 }
185
186 static float sb_fric_force_scale(Object *UNUSED(ob))
187 /* rescaling unit of drag [1 / sec] to somehow reasonable
188  * put it to a function here, so we can add user options later without touching simulation code
189  */
190 {
191         return (0.01f);
192 }
193
194 static float sb_time_scale(Object *ob)
195 /* defining the frames to *real* time relation */
196 {
197         SoftBody *sb= ob->soft; /* is supposed to be there */
198         if (sb) {
199                 return(sb->physics_speed);
200                 /* hrms .. this could be IPO as well :)
201                  * estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)]
202                  * 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions]  has period 65 frames
203                  * theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM)
204                  */
205         }
206         return (1.0f);
207         /*
208          * this would be frames/sec independent timing assuming 25 fps is default
209          * but does not work very well with NLA
210          * return (25.0f/scene->r.frs_sec)
211          */
212 }
213 /*--- frame based timing ---*/
214
215 /* helper functions for everything is animatable jow_go_for2_5 +++++++*/
216 /* introducing them here, because i know: steps in properties  ( at frame timing )
217  * will cause unwanted responses of the softbody system (which does inter frame calculations )
218  * so first 'cure' would be: interpolate linear in time ..
219  * Q: why do i write this?
220  * A: because it happened once, that some eger coder 'streamlined' code to fail.
221  * We DO linear interpolation for goals .. and i think we should do on animated properties as well
222  */
223
224 /* animate sb->maxgoal, sb->mingoal */
225 static float _final_goal(Object *ob, BodyPoint *bp)/*jow_go_for2_5 */
226 {
227         float f = -1999.99f;
228         if (ob) {
229                 SoftBody *sb= ob->soft; /* is supposed to be there */
230                 if (!(ob->softflag & OB_SB_GOAL)) return (0.0f);
231                 if (sb&&bp) {
232                         if (bp->goal < 0.0f) return (0.0f);
233                         f = sb->mingoal + bp->goal * fabsf(sb->maxgoal - sb->mingoal);
234                         f = pow(f, 4.0f);
235                         return (f);
236                 }
237         }
238         printf("_final_goal failed! sb or bp ==NULL\n");
239         return f; /*using crude but spot able values some times helps debuggin */
240 }
241
242 static float _final_mass(Object *ob, BodyPoint *bp)
243 {
244         if (ob) {
245                 SoftBody *sb= ob->soft; /* is supposed to be there */
246                 if (sb&&bp) {
247                         return(bp->mass*sb->nodemass);
248                 }
249         }
250         printf("_final_mass failed! sb or bp ==NULL\n");
251         return 1.0f;
252 }
253 /* helper functions for everything is animateble jow_go_for2_5 ------*/
254
255 /*+++ collider caching and dicing +++*/
256
257 /*
258  * for each target object/face the axis aligned bounding box (AABB) is stored
259  * faces parallel to global axes
260  * so only simple "value" in [min, max] ckecks are used
261  * float operations still
262  */
263
264 /* just an ID here to reduce the prob for killing objects
265  * ob->sumohandle points to we should not kill :)
266  */
267 static const int CCD_SAVETY = 190561;
268
269 typedef struct ccdf_minmax {
270         float minx, miny, minz, maxx, maxy, maxz;
271 } ccdf_minmax;
272
273
274 typedef struct ccd_Mesh {
275         int mvert_num, tri_num;
276         const MVert *mvert;
277         const MVert *mprevvert;
278         const MVertTri *tri;
279         int savety;
280         ccdf_minmax *mima;
281         /* Axis Aligned Bounding Box AABB */
282         float bbmin[3];
283         float bbmax[3];
284 } ccd_Mesh;
285
286
287 static ccd_Mesh *ccd_mesh_make(Object *ob)
288 {
289         CollisionModifierData *cmd;
290         ccd_Mesh *pccd_M = NULL;
291         ccdf_minmax *mima;
292         const MVertTri *vt;
293         float hull;
294         int i;
295
296         cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision);
297
298         /* first some paranoia checks */
299         if (!cmd) return NULL;
300         if (!cmd->mvert_num || !cmd->tri_num) return NULL;
301
302         pccd_M = MEM_mallocN(sizeof(ccd_Mesh), "ccd_Mesh");
303         pccd_M->mvert_num = cmd->mvert_num;
304         pccd_M->tri_num = cmd->tri_num;
305         pccd_M->savety  = CCD_SAVETY;
306         pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f;
307         pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f;
308         pccd_M->mprevvert=NULL;
309
310         /* blow it up with forcefield ranges */
311         hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft);
312
313         /* alloc and copy verts*/
314         pccd_M->mvert = MEM_dupallocN(cmd->xnew);
315         /* note that xnew coords are already in global space, */
316         /* determine the ortho BB */
317         for (i = 0; i < pccd_M->mvert_num; i++) {
318                 const float *v;
319
320                 /* evaluate limits */
321                 v = pccd_M->mvert[i].co;
322                 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
323                 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
324                 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
325
326                 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
327                 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
328                 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
329
330         }
331         /* alloc and copy faces*/
332         pccd_M->tri = MEM_dupallocN(cmd->tri);
333
334         /* OBBs for idea1 */
335         pccd_M->mima = MEM_mallocN(sizeof(ccdf_minmax) * pccd_M->tri_num, "ccd_Mesh_Faces_mima");
336
337
338         /* anyhoo we need to walk the list of faces and find OBB they live in */
339         for (i = 0, mima  = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
340                 const float *v;
341
342                 mima->minx=mima->miny=mima->minz=1e30f;
343                 mima->maxx=mima->maxy=mima->maxz=-1e30f;
344
345                 v = pccd_M->mvert[vt->tri[0]].co;
346                 mima->minx = min_ff(mima->minx, v[0] - hull);
347                 mima->miny = min_ff(mima->miny, v[1] - hull);
348                 mima->minz = min_ff(mima->minz, v[2] - hull);
349                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
350                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
351                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
352
353                 v = pccd_M->mvert[vt->tri[1]].co;
354                 mima->minx = min_ff(mima->minx, v[0] - hull);
355                 mima->miny = min_ff(mima->miny, v[1] - hull);
356                 mima->minz = min_ff(mima->minz, v[2] - hull);
357                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
358                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
359                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
360
361                 v = pccd_M->mvert[vt->tri[2]].co;
362                 mima->minx = min_ff(mima->minx, v[0] - hull);
363                 mima->miny = min_ff(mima->miny, v[1] - hull);
364                 mima->minz = min_ff(mima->minz, v[2] - hull);
365                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
366                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
367                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
368         }
369
370         return pccd_M;
371 }
372 static void ccd_mesh_update(Object *ob, ccd_Mesh *pccd_M)
373 {
374         CollisionModifierData *cmd;
375         ccdf_minmax *mima;
376         const MVertTri *vt;
377         float hull;
378         int i;
379
380         cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision);
381
382         /* first some paranoia checks */
383         if (!cmd) return;
384         if (!cmd->mvert_num || !cmd->tri_num) return;
385
386         if ((pccd_M->mvert_num != cmd->mvert_num) ||
387             (pccd_M->tri_num != cmd->tri_num))
388         {
389                 return;
390         }
391
392         pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f;
393         pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f;
394
395
396         /* blow it up with forcefield ranges */
397         hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft);
398
399         /* rotate current to previous */
400         if (pccd_M->mprevvert) MEM_freeN((void *)pccd_M->mprevvert);
401         pccd_M->mprevvert = pccd_M->mvert;
402         /* alloc and copy verts*/
403         pccd_M->mvert = MEM_dupallocN(cmd->xnew);
404         /* note that xnew coords are already in global space, */
405         /* determine the ortho BB */
406         for (i=0; i < pccd_M->mvert_num; i++) {
407                 const float *v;
408
409                 /* evaluate limits */
410                 v = pccd_M->mvert[i].co;
411                 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
412                 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
413                 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
414
415                 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
416                 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
417                 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
418
419                 /* evaluate limits */
420                 v = pccd_M->mprevvert[i].co;
421                 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
422                 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
423                 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
424
425                 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
426                 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
427                 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
428
429         }
430
431         /* anyhoo we need to walk the list of faces and find OBB they live in */
432         for (i = 0, mima  = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
433                 const float *v;
434
435                 mima->minx=mima->miny=mima->minz=1e30f;
436                 mima->maxx=mima->maxy=mima->maxz=-1e30f;
437
438                 /* mvert */
439                 v = pccd_M->mvert[vt->tri[0]].co;
440                 mima->minx = min_ff(mima->minx, v[0] - hull);
441                 mima->miny = min_ff(mima->miny, v[1] - hull);
442                 mima->minz = min_ff(mima->minz, v[2] - hull);
443                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
444                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
445                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
446
447                 v = pccd_M->mvert[vt->tri[1]].co;
448                 mima->minx = min_ff(mima->minx, v[0] - hull);
449                 mima->miny = min_ff(mima->miny, v[1] - hull);
450                 mima->minz = min_ff(mima->minz, v[2] - hull);
451                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
452                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
453                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
454
455                 v = pccd_M->mvert[vt->tri[2]].co;
456                 mima->minx = min_ff(mima->minx, v[0] - hull);
457                 mima->miny = min_ff(mima->miny, v[1] - hull);
458                 mima->minz = min_ff(mima->minz, v[2] - hull);
459                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
460                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
461                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
462
463
464                 /* mprevvert */
465                 v = pccd_M->mprevvert[vt->tri[0]].co;
466                 mima->minx = min_ff(mima->minx, v[0] - hull);
467                 mima->miny = min_ff(mima->miny, v[1] - hull);
468                 mima->minz = min_ff(mima->minz, v[2] - hull);
469                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
470                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
471                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
472
473                 v = pccd_M->mprevvert[vt->tri[1]].co;
474                 mima->minx = min_ff(mima->minx, v[0] - hull);
475                 mima->miny = min_ff(mima->miny, v[1] - hull);
476                 mima->minz = min_ff(mima->minz, v[2] - hull);
477                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
478                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
479                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
480
481                 v = pccd_M->mprevvert[vt->tri[2]].co;
482                 mima->minx = min_ff(mima->minx, v[0] - hull);
483                 mima->miny = min_ff(mima->miny, v[1] - hull);
484                 mima->minz = min_ff(mima->minz, v[2] - hull);
485                 mima->maxx = max_ff(mima->maxx, v[0] + hull);
486                 mima->maxy = max_ff(mima->maxy, v[1] + hull);
487                 mima->maxz = max_ff(mima->maxz, v[2] + hull);
488         }
489         return;
490 }
491
492 static void ccd_mesh_free(ccd_Mesh *ccdm)
493 {
494         if (ccdm && (ccdm->savety == CCD_SAVETY )) { /*make sure we're not nuking objects we don't know*/
495                 MEM_freeN((void *)ccdm->mvert);
496                 MEM_freeN((void *)ccdm->tri);
497                 if (ccdm->mprevvert) MEM_freeN((void *)ccdm->mprevvert);
498                 MEM_freeN(ccdm->mima);
499                 MEM_freeN(ccdm);
500                 ccdm = NULL;
501         }
502 }
503
504 static void ccd_build_deflector_hash_single(GHash *hash, Object *ob)
505 {
506         /* only with deflecting set */
507         if (ob->pd && ob->pd->deflect) {
508                 void **val_p;
509                 if (!BLI_ghash_ensure_p(hash, ob, &val_p)) {
510                         ccd_Mesh *ccdmesh = ccd_mesh_make(ob);
511                         *val_p = ccdmesh;
512                 }
513         }
514 }
515
516 /**
517  * \note collection overrides scene when not NULL.
518  */
519 static void ccd_build_deflector_hash(Depsgraph *depsgraph, Collection *collection, Object *vertexowner, GHash *hash)
520 {
521         if (!hash) return;
522
523         unsigned int numobjects;
524         Object **objects = BKE_collision_objects_create(depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
525
526         for (int i = 0; i < numobjects; i++) {
527                 Object *ob = objects[i];
528
529                 if (ob->type == OB_MESH) {
530                         ccd_build_deflector_hash_single(hash, ob);
531                 }
532         }
533
534         BKE_collision_objects_free(objects);
535 }
536
537 static void ccd_update_deflector_hash_single(GHash *hash, Object *ob)
538 {
539         if (ob->pd && ob->pd->deflect) {
540                 ccd_Mesh *ccdmesh = BLI_ghash_lookup(hash, ob);
541                 if (ccdmesh) {
542                         ccd_mesh_update(ob, ccdmesh);
543                 }
544         }
545 }
546
547 /**
548  * \note collection overrides scene when not NULL.
549  */
550 static void ccd_update_deflector_hash(Depsgraph *depsgraph, Collection *collection, Object *vertexowner, GHash *hash)
551 {
552         if ((!hash) || (!vertexowner)) return;
553
554         unsigned int numobjects;
555         Object **objects = BKE_collision_objects_create(depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
556
557         for (int i = 0; i < numobjects; i++) {
558                 Object *ob = objects[i];
559
560                 if (ob->type == OB_MESH) {
561                         ccd_update_deflector_hash_single(hash, ob);
562                 }
563         }
564
565         BKE_collision_objects_free(objects);
566 }
567
568 /*--- collider caching and dicing ---*/
569
570
571 static int count_mesh_quads(Mesh *me)
572 {
573         int a, result = 0;
574         const MPoly *mp = me->mpoly;
575
576         if (mp) {
577                 for (a = me->totpoly; a > 0; a--, mp++) {
578                         if (mp->totloop == 4) {
579                                 result++;
580                         }
581                 }
582         }
583         return result;
584 }
585
586 static void add_mesh_quad_diag_springs(Object *ob)
587 {
588         Mesh *me= ob->data;
589         /*BodyPoint *bp;*/ /*UNUSED*/
590         int a;
591
592         if (ob->soft) {
593                 int nofquads;
594                 //float s_shear = ob->soft->shearstiff*ob->soft->shearstiff;
595
596                 nofquads = count_mesh_quads(me);
597                 if (nofquads) {
598                         const MLoop *mloop = me->mloop;
599                         const MPoly *mp = me->mpoly;
600                         BodySpring *bs;
601
602                         /* resize spring-array to hold additional quad springs */
603                         ob->soft->bspring = MEM_recallocN(ob->soft->bspring, sizeof(BodySpring) * (ob->soft->totspring + nofquads * 2));
604
605                         /* fill the tail */
606                         a = 0;
607                         bs = &ob->soft->bspring[ob->soft->totspring];
608                         /*bp= ob->soft->bpoint; */ /*UNUSED*/
609                         for (a = me->totpoly; a > 0; a--, mp++) {
610                                 if (mp->totloop == 4) {
611                                         bs->v1 = mloop[mp->loopstart + 0].v;
612                                         bs->v2 = mloop[mp->loopstart + 2].v;
613                                         bs->springtype   = SB_STIFFQUAD;
614                                         bs++;
615                                         bs->v1 = mloop[mp->loopstart + 1].v;
616                                         bs->v2 = mloop[mp->loopstart + 3].v;
617                                         bs->springtype   = SB_STIFFQUAD;
618                                         bs++;
619                                 }
620                         }
621
622                         /* now we can announce new springs */
623                         ob->soft->totspring += nofquads * 2;
624                 }
625         }
626 }
627
628 static void add_2nd_order_roller(Object *ob, float UNUSED(stiffness), int *counter, int addsprings)
629 {
630         /*assume we have a softbody*/
631         SoftBody *sb= ob->soft; /* is supposed to be there */
632         BodyPoint *bp, *bpo;
633         BodySpring *bs, *bs2, *bs3= NULL;
634         int a, b, c, notthis= 0, v0;
635         if (!sb->bspring) {return;} /* we are 2nd order here so 1rst should have been build :) */
636         /* first run counting  second run adding */
637         *counter = 0;
638         if (addsprings) bs3 = ob->soft->bspring+ob->soft->totspring;
639         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
640                 /*scan for neighborhood*/
641                 bpo = NULL;
642                 v0  = (sb->totpoint-a);
643                 for (b=bp->nofsprings;b>0;b--) {
644                         bs = sb->bspring + bp->springs[b-1];
645                         /*nasty thing here that springs have two ends
646                         so here we have to make sure we examine the other */
647                         if (v0 == bs->v1) {
648                                 bpo = sb->bpoint+bs->v2;
649                                 notthis = bs->v2;
650                         }
651                         else {
652                                 if (v0 == bs->v2) {
653                                         bpo = sb->bpoint+bs->v1;
654                                         notthis = bs->v1;
655                                 }
656                                 else {
657                                         printf("oops we should not get here -  add_2nd_order_springs");
658                                 }
659                         }
660                         if (bpo) {/* so now we have a 2nd order humpdidump */
661                                 for (c=bpo->nofsprings;c>0;c--) {
662                                         bs2 = sb->bspring + bpo->springs[c-1];
663                                         if ((bs2->v1 != notthis) && (bs2->v1 > v0)) {
664                                                 (*counter)++;/*hit */
665                                                 if (addsprings) {
666                                                         bs3->v1= v0;
667                                                         bs3->v2= bs2->v1;
668                                                         bs3->springtype   = SB_BEND;
669                                                         bs3++;
670                                                 }
671                                         }
672                                         if ((bs2->v2 !=notthis)&&(bs2->v2 > v0)) {
673                                                 (*counter)++;  /* hit */
674                                                 if (addsprings) {
675                                                         bs3->v1= v0;
676                                                         bs3->v2= bs2->v2;
677                                                         bs3->springtype   = SB_BEND;
678                                                         bs3++;
679                                                 }
680
681                                         }
682                                 }
683
684                         }
685
686                 }
687                 /*scan for neighborhood done*/
688         }
689 }
690
691
692 static void add_2nd_order_springs(Object *ob, float stiffness)
693 {
694         int counter = 0;
695         BodySpring *bs_new;
696         stiffness *=stiffness;
697
698         add_2nd_order_roller(ob, stiffness, &counter, 0); /* counting */
699         if (counter) {
700                 /* resize spring-array to hold additional springs */
701                 bs_new= MEM_callocN((ob->soft->totspring + counter )*sizeof(BodySpring), "bodyspring");
702                 memcpy(bs_new, ob->soft->bspring, (ob->soft->totspring )*sizeof(BodySpring));
703
704                 if (ob->soft->bspring)
705                         MEM_freeN(ob->soft->bspring);
706                 ob->soft->bspring = bs_new;
707
708                 add_2nd_order_roller(ob, stiffness, &counter, 1); /* adding */
709                 ob->soft->totspring += counter;
710         }
711 }
712
713 static void add_bp_springlist(BodyPoint *bp, int springID)
714 {
715         int *newlist;
716
717         if (bp->springs == NULL) {
718                 bp->springs = MEM_callocN(sizeof(int), "bpsprings");
719                 bp->springs[0] = springID;
720                 bp->nofsprings = 1;
721         }
722         else {
723                 bp->nofsprings++;
724                 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings");
725                 memcpy(newlist, bp->springs, (bp->nofsprings-1)* sizeof(int));
726                 MEM_freeN(bp->springs);
727                 bp->springs = newlist;
728                 bp->springs[bp->nofsprings-1] = springID;
729         }
730 }
731
732 /* do this once when sb is build
733  * it is O(N^2) so scanning for springs every iteration is too expensive
734  */
735 static void build_bps_springlist(Object *ob)
736 {
737         SoftBody *sb= ob->soft; /* is supposed to be there */
738         BodyPoint *bp;
739         BodySpring *bs;
740         int a, b;
741
742         if (sb==NULL) return; /* paranoya check */
743
744         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
745                 /* throw away old list */
746                 if (bp->springs) {
747                         MEM_freeN(bp->springs);
748                         bp->springs=NULL;
749                 }
750                 /* scan for attached inner springs */
751                 for (b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) {
752                         if (( (sb->totpoint-a) == bs->v1) ) {
753                                 add_bp_springlist(bp, sb->totspring -b);
754                         }
755                         if (( (sb->totpoint-a) == bs->v2) ) {
756                                 add_bp_springlist(bp, sb->totspring -b);
757                         }
758                 }/*for springs*/
759         }/*for bp*/
760 }
761
762 static void calculate_collision_balls(Object *ob)
763 {
764         SoftBody *sb= ob->soft; /* is supposed to be there */
765         BodyPoint *bp;
766         BodySpring *bs;
767         int a, b, akku_count;
768         float min, max, akku;
769
770         if (sb==NULL) return; /* paranoya check */
771
772         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
773                 bp->colball=0;
774                 akku =0.0f;
775                 akku_count=0;
776                 min = 1e22f;
777                 max = -1e22f;
778                 /* first estimation based on attached */
779                 for (b=bp->nofsprings;b>0;b--) {
780                         bs = sb->bspring + bp->springs[b-1];
781                         if (bs->springtype == SB_EDGE) {
782                                 akku += bs->len;
783                                 akku_count++;
784                                 min = min_ff(bs->len, min);
785                                 max = max_ff(bs->len, max);
786                         }
787                 }
788
789                 if (akku_count > 0) {
790                         if (sb->sbc_mode == SBC_MODE_MANUAL) {
791                                 bp->colball=sb->colball;
792                         }
793                         if (sb->sbc_mode == SBC_MODE_AVG) {
794                                 bp->colball = akku/(float)akku_count*sb->colball;
795                         }
796                         if (sb->sbc_mode == SBC_MODE_MIN) {
797                                 bp->colball=min*sb->colball;
798                         }
799                         if (sb->sbc_mode == SBC_MODE_MAX) {
800                                 bp->colball=max*sb->colball;
801                         }
802                         if (sb->sbc_mode == SBC_MODE_AVGMINMAX) {
803                                 bp->colball = (min + max)/2.0f*sb->colball;
804                         }
805                 }
806                 else bp->colball=0;
807         }/*for bp*/
808 }
809
810
811 /* creates new softbody if didn't exist yet, makes new points and springs arrays */
812 static void renew_softbody(Scene *scene, Object *ob, int totpoint, int totspring)
813 {
814         SoftBody *sb;
815         int i;
816         short softflag;
817         if (ob->soft==NULL) ob->soft= sbNew(scene);
818         else free_softbody_intern(ob->soft);
819         sb= ob->soft;
820         softflag=ob->softflag;
821
822         if (totpoint) {
823                 sb->totpoint= totpoint;
824                 sb->totspring= totspring;
825
826                 sb->bpoint= MEM_mallocN(totpoint*sizeof(BodyPoint), "bodypoint");
827                 if (totspring)
828                         sb->bspring= MEM_mallocN(totspring*sizeof(BodySpring), "bodyspring");
829
830                         /* initialize BodyPoint array */
831                 for (i=0; i<totpoint; i++) {
832                         BodyPoint *bp = &sb->bpoint[i];
833
834
835                         /* hum as far as i see this is overridden by _final_goal() now jow_go_for2_5 */
836                         /* sadly breaks compatibility with older versions */
837                         /* but makes goals behave the same for meshes, lattices and curves */
838                         if (softflag & OB_SB_GOAL) {
839                                 bp->goal= sb->defgoal;
840                         }
841                         else {
842                                 bp->goal= 0.0f;
843                                 /* so this will definily be below SOFTGOALSNAP */
844                         }
845
846                         bp->nofsprings= 0;
847                         bp->springs= NULL;
848                         bp->choke = 0.0f;
849                         bp->choke2 = 0.0f;
850                         bp->frozen = 1.0f;
851                         bp->colball = 0.0f;
852                         bp->loc_flag = 0;
853                         bp->springweight = 1.0f;
854                         bp->mass = 1.0f;
855                 }
856         }
857 }
858
859 static void free_softbody_baked(SoftBody *sb)
860 {
861         SBVertex *key;
862         int k;
863
864         for (k=0; k<sb->totkey; k++) {
865                 key= *(sb->keys + k);
866                 if (key) MEM_freeN(key);
867         }
868         if (sb->keys) MEM_freeN(sb->keys);
869
870         sb->keys= NULL;
871         sb->totkey= 0;
872 }
873 static void free_scratch(SoftBody *sb)
874 {
875         if (sb->scratch) {
876                 /* todo make sure everything is cleaned up nicly */
877                 if (sb->scratch->colliderhash) {
878                         BLI_ghash_free(sb->scratch->colliderhash, NULL,
879                                         (GHashValFreeFP) ccd_mesh_free); /*this hoepfully will free all caches*/
880                         sb->scratch->colliderhash = NULL;
881                 }
882                 if (sb->scratch->bodyface) {
883                         MEM_freeN(sb->scratch->bodyface);
884                 }
885                 if (sb->scratch->Ref.ivert) {
886                         MEM_freeN(sb->scratch->Ref.ivert);
887                 }
888                 MEM_freeN(sb->scratch);
889                 sb->scratch = NULL;
890         }
891
892 }
893
894 /* only frees internal data */
895 static void free_softbody_intern(SoftBody *sb)
896 {
897         if (sb) {
898                 int a;
899                 BodyPoint *bp;
900
901                 if (sb->bpoint) {
902                         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
903                                 /* free spring list */
904                                 if (bp->springs != NULL) {
905                                         MEM_freeN(bp->springs);
906                                 }
907                         }
908                         MEM_freeN(sb->bpoint);
909                 }
910
911                 if (sb->bspring) MEM_freeN(sb->bspring);
912
913                 sb->totpoint= sb->totspring= 0;
914                 sb->bpoint= NULL;
915                 sb->bspring= NULL;
916
917                 free_scratch(sb);
918                 free_softbody_baked(sb);
919         }
920 }
921
922
923 /* ************ dynamics ********** */
924
925 /* the most general (micro physics correct) way to do collision
926  * (only needs the current particle position)
927  *
928  * it actually checks if the particle intrudes a short range force field generated
929  * by the faces of the target object and returns a force to drive the particel out
930  * the strength of the field grows exponetially if the particle is on the 'wrong' side of the face
931  * 'wrong' side : projection to the face normal is negative (all referred to a vertex in the face)
932  *
933  * flaw of this: 'fast' particles as well as 'fast' colliding faces
934  * give a 'tunnel' effect such that the particle passes through the force field
935  * without ever 'seeing' it
936  * this is fully compliant to heisenberg: h >= fuzzy(location) * fuzzy(time)
937  * besides our h is way larger than in QM because forces propagate way slower here
938  * we have to deal with fuzzy(time) in the range of 1/25 seconds (typical frame rate)
939  * yup collision targets are not known here any better
940  * and 1/25 second is looong compared to real collision events
941  * Q: why not use 'simple' collision here like bouncing back a particle
942  *   --> reverting is velocity on the face normal
943  * A: because our particles are not alone here
944  *    and need to tell their neighbors exactly what happens via spring forces
945  * unless sbObjectStep( .. ) is called on sub frame timing level
946  * BTW that also questions the use of a 'implicit' solvers on softbodies
947  * since that would only valid for 'slow' moving collision targets and dito particles
948  */
949
950 /* +++ dependency information functions*/
951
952 /**
953  * \note collection overrides scene when not NULL.
954  */
955 static int query_external_colliders(Depsgraph *depsgraph, Collection *collection)
956 {
957         unsigned int numobjects;
958         Object **objects = BKE_collision_objects_create(depsgraph, NULL, collection, &numobjects, eModifierType_Collision);
959         BKE_collision_objects_free(objects);
960
961         return (numobjects != 0);
962 }
963 /* --- dependency information functions*/
964
965
966 /* +++ the aabb "force" section*/
967 static int sb_detect_aabb_collisionCached(float UNUSED(force[3]), struct Object *vertexowner, float UNUSED(time))
968 {
969         Object *ob;
970         SoftBody *sb=vertexowner->soft;
971         GHash *hash;
972         GHashIterator *ihash;
973         float  aabbmin[3], aabbmax[3];
974         int deflected=0;
975
976         if ((sb == NULL) || (sb->scratch ==NULL)) return 0;
977         copy_v3_v3(aabbmin, sb->scratch->aabbmin);
978         copy_v3_v3(aabbmax, sb->scratch->aabbmax);
979
980         hash  = vertexowner->soft->scratch->colliderhash;
981         ihash = BLI_ghashIterator_new(hash);
982         while (!BLI_ghashIterator_done(ihash)) {
983
984                 ccd_Mesh *ccdm = BLI_ghashIterator_getValue     (ihash);
985                 ob             = BLI_ghashIterator_getKey       (ihash);
986                 {
987                         /* only with deflecting set */
988                         if (ob->pd && ob->pd->deflect) {
989                                 if (ccdm) {
990                                         if ((aabbmax[0] < ccdm->bbmin[0]) ||
991                                             (aabbmax[1] < ccdm->bbmin[1]) ||
992                                             (aabbmax[2] < ccdm->bbmin[2]) ||
993                                             (aabbmin[0] > ccdm->bbmax[0]) ||
994                                             (aabbmin[1] > ccdm->bbmax[1]) ||
995                                             (aabbmin[2] > ccdm->bbmax[2]) )
996                                         {
997                                                 /* boxes don't intersect */
998                                                 BLI_ghashIterator_step(ihash);
999                                                 continue;
1000                                         }
1001
1002                                         /* so now we have the 2 boxes overlapping */
1003                                         /* forces actually not used */
1004                                         deflected = 2;
1005
1006                                 }
1007                                 else {
1008                                         /*aye that should be cached*/
1009                                         printf("missing cache error\n");
1010                                         BLI_ghashIterator_step(ihash);
1011                                         continue;
1012                                 }
1013                         } /* if (ob->pd && ob->pd->deflect) */
1014                         BLI_ghashIterator_step(ihash);
1015                 }
1016         } /* while () */
1017         BLI_ghashIterator_free(ihash);
1018         return deflected;
1019 }
1020 /* --- the aabb section*/
1021
1022
1023 /* +++ the face external section*/
1024 static int sb_detect_face_pointCached(float face_v1[3], float face_v2[3], float face_v3[3], float *damp,
1025                                       float force[3], struct Object *vertexowner, float time)
1026 {
1027         Object *ob;
1028         GHash *hash;
1029         GHashIterator *ihash;
1030         float nv1[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1031         float facedist, outerfacethickness, tune = 10.f;
1032         int a, deflected=0;
1033
1034         aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]);
1035         aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]);
1036         aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]);
1037         aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]);
1038         aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]);
1039         aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]);
1040
1041         /* calculate face normal once again SIGH */
1042         sub_v3_v3v3(edge1, face_v1, face_v2);
1043         sub_v3_v3v3(edge2, face_v3, face_v2);
1044         cross_v3_v3v3(d_nvect, edge2, edge1);
1045         normalize_v3(d_nvect);
1046
1047
1048         hash  = vertexowner->soft->scratch->colliderhash;
1049         ihash = BLI_ghashIterator_new(hash);
1050         while (!BLI_ghashIterator_done(ihash)) {
1051
1052                 ccd_Mesh *ccdm = BLI_ghashIterator_getValue     (ihash);
1053                 ob             = BLI_ghashIterator_getKey       (ihash);
1054                 {
1055                         /* only with deflecting set */
1056                         if (ob->pd && ob->pd->deflect) {
1057                                 const MVert *mvert= NULL;
1058                                 const MVert *mprevvert= NULL;
1059                                 if (ccdm) {
1060                                         mvert = ccdm->mvert;
1061                                         a     = ccdm->mvert_num;
1062                                         mprevvert= ccdm->mprevvert;
1063                                         outerfacethickness = ob->pd->pdef_sboft;
1064                                         if ((aabbmax[0] < ccdm->bbmin[0]) ||
1065                                             (aabbmax[1] < ccdm->bbmin[1]) ||
1066                                             (aabbmax[2] < ccdm->bbmin[2]) ||
1067                                             (aabbmin[0] > ccdm->bbmax[0]) ||
1068                                             (aabbmin[1] > ccdm->bbmax[1]) ||
1069                                             (aabbmin[2] > ccdm->bbmax[2]) )
1070                                         {
1071                                                 /* boxes don't intersect */
1072                                                 BLI_ghashIterator_step(ihash);
1073                                                 continue;
1074                                         }
1075
1076                                 }
1077                                 else {
1078                                         /*aye that should be cached*/
1079                                         printf("missing cache error\n");
1080                                         BLI_ghashIterator_step(ihash);
1081                                         continue;
1082                                 }
1083
1084
1085                                 /* use mesh*/
1086                                 if (mvert) {
1087                                         while (a) {
1088                                                 copy_v3_v3(nv1, mvert[a-1].co);
1089                                                 if (mprevvert) {
1090                                                         mul_v3_fl(nv1, time);
1091                                                         madd_v3_v3fl(nv1, mprevvert[a - 1].co, 1.0f - time);
1092                                                 }
1093                                                 /* origin to face_v2*/
1094                                                 sub_v3_v3(nv1, face_v2);
1095                                                 facedist = dot_v3v3(nv1, d_nvect);
1096                                                 if (ABS(facedist)<outerfacethickness) {
1097                                                         if (isect_point_tri_prism_v3(nv1, face_v1, face_v2, face_v3) ) {
1098                                                                 float df;
1099                                                                 if (facedist > 0) {
1100                                                                         df = (outerfacethickness-facedist)/outerfacethickness;
1101                                                                 }
1102                                                                 else {
1103                                                                         df = (outerfacethickness+facedist)/outerfacethickness;
1104                                                                 }
1105
1106                                                                 *damp=df*tune*ob->pd->pdef_sbdamp;
1107
1108                                                                 df = 0.01f * expf(-100.0f * df);
1109                                                                 madd_v3_v3fl(force, d_nvect, -df);
1110                                                                 deflected = 3;
1111                                                         }
1112                                                 }
1113                                                 a--;
1114                                         }/* while (a)*/
1115                                 } /* if (mvert) */
1116                         } /* if (ob->pd && ob->pd->deflect) */
1117                         BLI_ghashIterator_step(ihash);
1118                 }
1119         } /* while () */
1120         BLI_ghashIterator_free(ihash);
1121         return deflected;
1122 }
1123
1124
1125 static int sb_detect_face_collisionCached(float face_v1[3], float face_v2[3], float face_v3[3], float *damp,
1126                                           float force[3], struct Object *vertexowner, float time)
1127 {
1128         Object *ob;
1129         GHash *hash;
1130         GHashIterator *ihash;
1131         float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1132         float t, tune = 10.0f;
1133         int a, deflected=0;
1134
1135         aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]);
1136         aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]);
1137         aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]);
1138         aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]);
1139         aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]);
1140         aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]);
1141
1142         hash  = vertexowner->soft->scratch->colliderhash;
1143         ihash = BLI_ghashIterator_new(hash);
1144         while (!BLI_ghashIterator_done(ihash)) {
1145
1146                 ccd_Mesh *ccdm = BLI_ghashIterator_getValue     (ihash);
1147                 ob             = BLI_ghashIterator_getKey       (ihash);
1148                 {
1149                         /* only with deflecting set */
1150                         if (ob->pd && ob->pd->deflect) {
1151                                 const MVert *mvert = NULL;
1152                                 const MVert *mprevvert = NULL;
1153                                 const MVertTri *vt = NULL;
1154                                 const ccdf_minmax *mima = NULL;
1155
1156                                 if (ccdm) {
1157                                         mvert = ccdm->mvert;
1158                                         vt = ccdm->tri;
1159                                         mprevvert = ccdm->mprevvert;
1160                                         mima = ccdm->mima;
1161                                         a = ccdm->tri_num;
1162
1163                                         if ((aabbmax[0] < ccdm->bbmin[0]) ||
1164                                             (aabbmax[1] < ccdm->bbmin[1]) ||
1165                                             (aabbmax[2] < ccdm->bbmin[2]) ||
1166                                             (aabbmin[0] > ccdm->bbmax[0]) ||
1167                                             (aabbmin[1] > ccdm->bbmax[1]) ||
1168                                             (aabbmin[2] > ccdm->bbmax[2]) )
1169                                         {
1170                                                 /* boxes don't intersect */
1171                                                 BLI_ghashIterator_step(ihash);
1172                                                 continue;
1173                                         }
1174
1175                                 }
1176                                 else {
1177                                         /*aye that should be cached*/
1178                                         printf("missing cache error\n");
1179                                         BLI_ghashIterator_step(ihash);
1180                                         continue;
1181                                 }
1182
1183
1184                                 /* use mesh*/
1185                                 while (a--) {
1186                                         if ((aabbmax[0] < mima->minx) ||
1187                                             (aabbmin[0] > mima->maxx) ||
1188                                             (aabbmax[1] < mima->miny) ||
1189                                             (aabbmin[1] > mima->maxy) ||
1190                                             (aabbmax[2] < mima->minz) ||
1191                                             (aabbmin[2] > mima->maxz))
1192                                         {
1193                                                 mima++;
1194                                                 vt++;
1195                                                 continue;
1196                                         }
1197
1198
1199                                         if (mvert) {
1200
1201                                                 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1202                                                 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1203                                                 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1204
1205                                                 if (mprevvert) {
1206                                                         mul_v3_fl(nv1, time);
1207                                                         madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1208
1209                                                         mul_v3_fl(nv2, time);
1210                                                         madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1211
1212                                                         mul_v3_fl(nv3, time);
1213                                                         madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1214                                                 }
1215                                         }
1216
1217                                         /* switch origin to be nv2*/
1218                                         sub_v3_v3v3(edge1, nv1, nv2);
1219                                         sub_v3_v3v3(edge2, nv3, nv2);
1220                                         cross_v3_v3v3(d_nvect, edge2, edge1);
1221                                         normalize_v3(d_nvect);
1222                                         if (isect_line_segment_tri_v3(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) ||
1223                                             isect_line_segment_tri_v3(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) ||
1224                                             isect_line_segment_tri_v3(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL) )
1225                                         {
1226                                                 madd_v3_v3fl(force, d_nvect, -0.5f);
1227                                                 *damp=tune*ob->pd->pdef_sbdamp;
1228                                                 deflected = 2;
1229                                         }
1230                                         mima++;
1231                                         vt++;
1232                                 }/* while a */
1233                         } /* if (ob->pd && ob->pd->deflect) */
1234                         BLI_ghashIterator_step(ihash);
1235                 }
1236         } /* while () */
1237         BLI_ghashIterator_free(ihash);
1238         return deflected;
1239 }
1240
1241
1242
1243 static void scan_for_ext_face_forces(Object *ob, float timenow)
1244 {
1245         SoftBody *sb = ob->soft;
1246         BodyFace *bf;
1247         int a;
1248         float damp=0.0f, choke=1.0f;
1249         float tune = -10.0f;
1250         float feedback[3];
1251
1252         if (sb && sb->scratch->totface) {
1253
1254
1255                 bf = sb->scratch->bodyface;
1256                 for (a=0; a<sb->scratch->totface; a++, bf++) {
1257                         bf->ext_force[0]=bf->ext_force[1]=bf->ext_force[2]=0.0f;
1258 /*+++edges intruding*/
1259                         bf->flag &= ~BFF_INTERSECT;
1260                         zero_v3(feedback);
1261                         if (sb_detect_face_collisionCached(
1262                                 sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos,
1263                                 &damp, feedback, ob, timenow))
1264                         {
1265                                 madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
1266                                 madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
1267                                 madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
1268 //                              madd_v3_v3fl(bf->ext_force, feedback, tune);
1269                                 bf->flag |= BFF_INTERSECT;
1270                                 choke = min_ff(max_ff(damp, choke), 1.0f);
1271                         }
1272 /*---edges intruding*/
1273
1274 /*+++ close vertices*/
1275                         if (( bf->flag & BFF_INTERSECT)==0) {
1276                                 bf->flag &= ~BFF_CLOSEVERT;
1277                                 tune = -1.0f;
1278                                 zero_v3(feedback);
1279                                 if (sb_detect_face_pointCached(
1280                                         sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos,
1281                                         &damp,  feedback, ob, timenow))
1282                                 {
1283                                         madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
1284                                         madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
1285                                         madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
1286 //                                      madd_v3_v3fl(bf->ext_force, feedback, tune);
1287                                         bf->flag |= BFF_CLOSEVERT;
1288                                         choke = min_ff(max_ff(damp, choke), 1.0f);
1289                                 }
1290                         }
1291 /*--- close vertices*/
1292                 }
1293                 bf = sb->scratch->bodyface;
1294                 for (a=0; a<sb->scratch->totface; a++, bf++) {
1295                         if (( bf->flag & BFF_INTERSECT) || ( bf->flag & BFF_CLOSEVERT)) {
1296                                 sb->bpoint[bf->v1].choke2 = max_ff(sb->bpoint[bf->v1].choke2, choke);
1297                                 sb->bpoint[bf->v2].choke2 = max_ff(sb->bpoint[bf->v2].choke2, choke);
1298                                 sb->bpoint[bf->v3].choke2 = max_ff(sb->bpoint[bf->v3].choke2, choke);
1299                         }
1300                 }
1301         }
1302 }
1303
1304 /*  --- the face external section*/
1305
1306
1307 /* +++ the spring external section*/
1308
1309 static int sb_detect_edge_collisionCached(float edge_v1[3], float edge_v2[3], float *damp,
1310                                                                    float force[3], struct Object *vertexowner, float time)
1311 {
1312         Object *ob;
1313         GHash *hash;
1314         GHashIterator *ihash;
1315         float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1316         float t, el;
1317         int a, deflected=0;
1318
1319         minmax_v3v3_v3(aabbmin, aabbmax, edge_v1);
1320         minmax_v3v3_v3(aabbmin, aabbmax, edge_v2);
1321
1322         el = len_v3v3(edge_v1, edge_v2);
1323
1324         hash  = vertexowner->soft->scratch->colliderhash;
1325         ihash = BLI_ghashIterator_new(hash);
1326         while (!BLI_ghashIterator_done(ihash)) {
1327
1328                 ccd_Mesh *ccdm = BLI_ghashIterator_getValue     (ihash);
1329                 ob             = BLI_ghashIterator_getKey       (ihash);
1330                 {
1331                         /* only with deflecting set */
1332                         if (ob->pd && ob->pd->deflect) {
1333                                 const MVert *mvert = NULL;
1334                                 const MVert *mprevvert = NULL;
1335                                 const MVertTri *vt = NULL;
1336                                 const ccdf_minmax *mima = NULL;
1337
1338                                 if (ccdm) {
1339                                         mvert = ccdm->mvert;
1340                                         mprevvert = ccdm->mprevvert;
1341                                         vt = ccdm->tri;
1342                                         mima = ccdm->mima;
1343                                         a = ccdm->tri_num;
1344
1345                                         if ((aabbmax[0] < ccdm->bbmin[0]) ||
1346                                             (aabbmax[1] < ccdm->bbmin[1]) ||
1347                                             (aabbmax[2] < ccdm->bbmin[2]) ||
1348                                             (aabbmin[0] > ccdm->bbmax[0]) ||
1349                                             (aabbmin[1] > ccdm->bbmax[1]) ||
1350                                             (aabbmin[2] > ccdm->bbmax[2]) )
1351                                         {
1352                                                 /* boxes don't intersect */
1353                                                 BLI_ghashIterator_step(ihash);
1354                                                 continue;
1355                                         }
1356
1357                                 }
1358                                 else {
1359                                         /*aye that should be cached*/
1360                                         printf("missing cache error\n");
1361                                         BLI_ghashIterator_step(ihash);
1362                                         continue;
1363                                 }
1364
1365
1366                                 /* use mesh*/
1367                                 while (a--) {
1368                                         if ((aabbmax[0] < mima->minx) ||
1369                                             (aabbmin[0] > mima->maxx) ||
1370                                             (aabbmax[1] < mima->miny) ||
1371                                             (aabbmin[1] > mima->maxy) ||
1372                                             (aabbmax[2] < mima->minz) ||
1373                                             (aabbmin[2] > mima->maxz))
1374                                         {
1375                                                 mima++;
1376                                                 vt++;
1377                                                 continue;
1378                                         }
1379
1380
1381                                         if (mvert) {
1382
1383                                                 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1384                                                 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1385                                                 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1386
1387                                                 if (mprevvert) {
1388                                                         mul_v3_fl(nv1, time);
1389                                                         madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1390
1391                                                         mul_v3_fl(nv2, time);
1392                                                         madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1393
1394                                                         mul_v3_fl(nv3, time);
1395                                                         madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1396                                                 }
1397                                         }
1398
1399                                         /* switch origin to be nv2*/
1400                                         sub_v3_v3v3(edge1, nv1, nv2);
1401                                         sub_v3_v3v3(edge2, nv3, nv2);
1402
1403                                         cross_v3_v3v3(d_nvect, edge2, edge1);
1404                                         normalize_v3(d_nvect);
1405                                         if (isect_line_segment_tri_v3(edge_v1, edge_v2, nv1, nv2, nv3, &t, NULL)) {
1406                                                 float v1[3], v2[3];
1407                                                 float intrusiondepth, i1, i2;
1408                                                 sub_v3_v3v3(v1, edge_v1, nv2);
1409                                                 sub_v3_v3v3(v2, edge_v2, nv2);
1410                                                 i1 = dot_v3v3(v1, d_nvect);
1411                                                 i2 = dot_v3v3(v2, d_nvect);
1412                                                 intrusiondepth = -min_ff(i1, i2) / el;
1413                                                 madd_v3_v3fl(force, d_nvect, intrusiondepth);
1414                                                 *damp=ob->pd->pdef_sbdamp;
1415                                                 deflected = 2;
1416                                         }
1417
1418                                         mima++;
1419                                         vt++;
1420                                 }/* while a */
1421                         } /* if (ob->pd && ob->pd->deflect) */
1422                         BLI_ghashIterator_step(ihash);
1423                 }
1424         } /* while () */
1425         BLI_ghashIterator_free(ihash);
1426         return deflected;
1427 }
1428
1429 static void _scan_for_ext_spring_forces(Scene *scene, Object *ob, float timenow, int ifirst, int ilast, struct ListBase *effectors)
1430 {
1431         SoftBody *sb = ob->soft;
1432         int a;
1433         float damp;
1434         float feedback[3];
1435
1436         if (sb && sb->totspring) {
1437                 for (a=ifirst; a<ilast; a++) {
1438                         BodySpring *bs = &sb->bspring[a];
1439                         bs->ext_force[0]=bs->ext_force[1]=bs->ext_force[2]=0.0f;
1440                         feedback[0]=feedback[1]=feedback[2]=0.0f;
1441                         bs->flag &= ~BSF_INTERSECT;
1442
1443                         if (bs->springtype == SB_EDGE) {
1444                                 /* +++ springs colliding */
1445                                 if (ob->softflag & OB_SB_EDGECOLL) {
1446                                         if ( sb_detect_edge_collisionCached (sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos,
1447                                                 &damp, feedback, ob, timenow)) {
1448                                                         add_v3_v3(bs->ext_force, feedback);
1449                                                         bs->flag |= BSF_INTERSECT;
1450                                                         //bs->cf=damp;
1451                                                         bs->cf=sb->choke*0.01f;
1452
1453                                         }
1454                                 }
1455                                 /* ---- springs colliding */
1456
1457                                 /* +++ springs seeing wind ... n stuff depending on their orientation*/
1458                                 /* note we don't use sb->mediafrict but use sb->aeroedge for magnitude of effect*/
1459                                 if (sb->aeroedge) {
1460                                         float vel[3], sp[3], pr[3], force[3];
1461                                         float f, windfactor  = 0.25f;
1462                                         /*see if we have wind*/
1463                                         if (effectors) {
1464                                                 EffectedPoint epoint;
1465                                                 float speed[3] = {0.0f, 0.0f, 0.0f};
1466                                                 float pos[3];
1467                                                 mid_v3_v3v3(pos, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos);
1468                                                 mid_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec);
1469                                                 pd_point_from_soft(scene, pos, vel, -1, &epoint);
1470                                                 BKE_effectors_apply(effectors, NULL, sb->effector_weights, &epoint, force, speed);
1471
1472                                                 mul_v3_fl(speed, windfactor);
1473                                                 add_v3_v3(vel, speed);
1474                                         }
1475                                         /* media in rest */
1476                                         else {
1477                                                 add_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec);
1478                                         }
1479                                         f = normalize_v3(vel);
1480                                         f = -0.0001f*f*f*sb->aeroedge;
1481                                         /* (todo) add a nice angle dependent function done for now BUT */
1482                                         /* still there could be some nice drag/lift function, but who needs it */
1483
1484                                         sub_v3_v3v3(sp, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos);
1485                                         project_v3_v3v3(pr, vel, sp);
1486                                         sub_v3_v3(vel, pr);
1487                                         normalize_v3(vel);
1488                                         if (ob->softflag & OB_SB_AERO_ANGLE) {
1489                                                 normalize_v3(sp);
1490                                                 madd_v3_v3fl(bs->ext_force, vel, f * (1.0f - fabsf(dot_v3v3(vel, sp))));
1491                                         }
1492                                         else {
1493                                                 madd_v3_v3fl(bs->ext_force, vel, f); // to keep compatible with 2.45 release files
1494                                         }
1495                                 }
1496                                 /* --- springs seeing wind */
1497                         }
1498                 }
1499         }
1500 }
1501
1502
1503 static void scan_for_ext_spring_forces(struct Depsgraph *depsgraph, Scene *scene, Object *ob, float timenow)
1504 {
1505         SoftBody *sb = ob->soft;
1506
1507         ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, sb->effector_weights);
1508         _scan_for_ext_spring_forces(scene, ob, timenow, 0, sb->totspring, effectors);
1509         BKE_effectors_free(effectors);
1510 }
1511
1512 static void *exec_scan_for_ext_spring_forces(void *data)
1513 {
1514         SB_thread_context *pctx = (SB_thread_context*)data;
1515         _scan_for_ext_spring_forces(pctx->scene, pctx->ob, pctx->timenow, pctx->ifirst, pctx->ilast, pctx->effectors);
1516         return NULL;
1517 }
1518
1519 static void sb_sfesf_threads_run(struct Depsgraph *depsgraph, Scene *scene, struct Object *ob, float timenow, int totsprings, int *UNUSED(ptr_to_break_func(void)))
1520 {
1521         ListBase threads;
1522         SB_thread_context *sb_threads;
1523         int i, totthread, left, dec;
1524         int lowsprings =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */
1525
1526         ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, ob->soft->effector_weights);
1527
1528         /* figure the number of threads while preventing pretty pointless threading overhead */
1529         totthread= BKE_scene_num_threads(scene);
1530         /* what if we got zillions of CPUs running but less to spread*/
1531         while ((totsprings/totthread < lowsprings) && (totthread > 1)) {
1532                 totthread--;
1533         }
1534
1535         sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBSpringsThread");
1536         memset(sb_threads, 0, sizeof(SB_thread_context)*totthread);
1537         left = totsprings;
1538         dec = totsprings/totthread +1;
1539         for (i=0; i<totthread; i++) {
1540                 sb_threads[i].scene = scene;
1541                 sb_threads[i].ob = ob;
1542                 sb_threads[i].forcetime = 0.0; // not used here
1543                 sb_threads[i].timenow = timenow;
1544                 sb_threads[i].ilast   = left;
1545                 left = left - dec;
1546                 if (left >0) {
1547                         sb_threads[i].ifirst  = left;
1548                 }
1549                 else
1550                         sb_threads[i].ifirst  = 0;
1551                 sb_threads[i].effectors = effectors;
1552                 sb_threads[i].do_deflector = false;// not used here
1553                 sb_threads[i].fieldfactor = 0.0f;// not used here
1554                 sb_threads[i].windfactor  = 0.0f;// not used here
1555                 sb_threads[i].nr= i;
1556                 sb_threads[i].tot= totthread;
1557         }
1558         if (totthread > 1) {
1559                 BLI_threadpool_init(&threads, exec_scan_for_ext_spring_forces, totthread);
1560
1561                 for (i=0; i<totthread; i++)
1562                         BLI_threadpool_insert(&threads, &sb_threads[i]);
1563
1564                 BLI_threadpool_end(&threads);
1565         }
1566         else
1567                 exec_scan_for_ext_spring_forces(&sb_threads[0]);
1568         /* clean up */
1569         MEM_freeN(sb_threads);
1570
1571         BKE_effectors_free(effectors);
1572 }
1573
1574
1575 /* --- the spring external section*/
1576
1577 static int choose_winner(float*w, float* pos, float*a, float*b, float*c, float*ca, float*cb, float*cc)
1578 {
1579         float mindist, cp;
1580         int winner =1;
1581         mindist = fabsf(dot_v3v3(pos, a));
1582
1583         cp = fabsf(dot_v3v3(pos, b));
1584         if ( mindist < cp ) {
1585                 mindist = cp;
1586                 winner =2;
1587         }
1588
1589         cp = fabsf(dot_v3v3(pos, c));
1590         if (mindist < cp ) {
1591                 mindist = cp;
1592                 winner =3;
1593         }
1594         switch (winner) {
1595                 case 1: copy_v3_v3(w, ca); break;
1596                 case 2: copy_v3_v3(w, cb); break;
1597                 case 3: copy_v3_v3(w, cc);
1598         }
1599         return(winner);
1600 }
1601
1602
1603
1604 static int sb_detect_vertex_collisionCached(
1605         float opco[3], float facenormal[3], float *damp,
1606         float force[3], struct Object *vertexowner,
1607         float time, float vel[3], float *intrusion)
1608 {
1609         Object *ob= NULL;
1610         GHash *hash;
1611         GHashIterator *ihash;
1612         float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], dv1[3], ve[3], avel[3] = {0.0, 0.0, 0.0},
1613               vv1[3], vv2[3], vv3[3], coledge[3] = {0.0f, 0.0f, 0.0f}, mindistedge = 1000.0f,
1614               outerforceaccu[3], innerforceaccu[3],
1615               facedist, /* n_mag, */ /* UNUSED */ force_mag_norm, minx, miny, minz, maxx, maxy, maxz,
1616               innerfacethickness = -0.5f, outerfacethickness = 0.2f,
1617               ee = 5.0f, ff = 0.1f, fa=1;
1618         int a, deflected=0, cavel=0, ci=0;
1619 /* init */
1620         *intrusion = 0.0f;
1621         hash  = vertexowner->soft->scratch->colliderhash;
1622         ihash = BLI_ghashIterator_new(hash);
1623         outerforceaccu[0]=outerforceaccu[1]=outerforceaccu[2]=0.0f;
1624         innerforceaccu[0]=innerforceaccu[1]=innerforceaccu[2]=0.0f;
1625 /* go */
1626         while (!BLI_ghashIterator_done(ihash)) {
1627
1628                 ccd_Mesh *ccdm = BLI_ghashIterator_getValue     (ihash);
1629                 ob             = BLI_ghashIterator_getKey       (ihash);
1630                 {
1631                         /* only with deflecting set */
1632                         if (ob->pd && ob->pd->deflect) {
1633                                 const MVert *mvert = NULL;
1634                                 const MVert *mprevvert = NULL;
1635                                 const MVertTri *vt = NULL;
1636                                 const ccdf_minmax *mima = NULL;
1637
1638                                 if (ccdm) {
1639                                         mvert = ccdm->mvert;
1640                                         mprevvert = ccdm->mprevvert;
1641                                         vt = ccdm->tri;
1642                                         mima = ccdm->mima;
1643                                         a = ccdm->tri_num;
1644
1645                                         minx = ccdm->bbmin[0];
1646                                         miny = ccdm->bbmin[1];
1647                                         minz = ccdm->bbmin[2];
1648
1649                                         maxx = ccdm->bbmax[0];
1650                                         maxy = ccdm->bbmax[1];
1651                                         maxz = ccdm->bbmax[2];
1652
1653                                         if ((opco[0] < minx) ||
1654                                             (opco[1] < miny) ||
1655                                             (opco[2] < minz) ||
1656                                             (opco[0] > maxx) ||
1657                                             (opco[1] > maxy) ||
1658                                             (opco[2] > maxz) )
1659                                         {
1660                                                 /* outside the padded boundbox --> collision object is too far away */
1661                                                 BLI_ghashIterator_step(ihash);
1662                                                 continue;
1663                                         }
1664                                 }
1665                                 else {
1666                                         /*aye that should be cached*/
1667                                         printf("missing cache error\n");
1668                                         BLI_ghashIterator_step(ihash);
1669                                         continue;
1670                                 }
1671
1672                                 /* do object level stuff */
1673                                 /* need to have user control for that since it depends on model scale */
1674                                 innerfacethickness = -ob->pd->pdef_sbift;
1675                                 outerfacethickness =  ob->pd->pdef_sboft;
1676                                 fa = (ff*outerfacethickness-outerfacethickness);
1677                                 fa *= fa;
1678                                 fa = 1.0f/fa;
1679                                 avel[0]=avel[1]=avel[2]=0.0f;
1680                                 /* use mesh*/
1681                                 while (a--) {
1682                                         if ((opco[0] < mima->minx) ||
1683                                             (opco[0] > mima->maxx) ||
1684                                             (opco[1] < mima->miny) ||
1685                                             (opco[1] > mima->maxy) ||
1686                                             (opco[2] < mima->minz) ||
1687                                             (opco[2] > mima->maxz))
1688                                         {
1689                                                 mima++;
1690                                                 vt++;
1691                                                 continue;
1692                                         }
1693
1694                                         if (mvert) {
1695
1696                                                 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1697                                                 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1698                                                 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1699
1700                                                 if (mprevvert) {
1701                                                         /* grab the average speed of the collider vertices
1702                                                          * before we spoil nvX
1703                                                          * humm could be done once a SB steps but then we' need to store that too
1704                                                          * since the AABB reduced propabitlty to get here drasticallly
1705                                                          * it might be a nice tradeof CPU <--> memory
1706                                                          */
1707                                                         sub_v3_v3v3(vv1, nv1, mprevvert[vt->tri[0]].co);
1708                                                         sub_v3_v3v3(vv2, nv2, mprevvert[vt->tri[1]].co);
1709                                                         sub_v3_v3v3(vv3, nv3, mprevvert[vt->tri[2]].co);
1710
1711                                                         mul_v3_fl(nv1, time);
1712                                                         madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1713
1714                                                         mul_v3_fl(nv2, time);
1715                                                         madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1716
1717                                                         mul_v3_fl(nv3, time);
1718                                                         madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1719                                                 }
1720                                         }
1721
1722                                         /* switch origin to be nv2*/
1723                                         sub_v3_v3v3(edge1, nv1, nv2);
1724                                         sub_v3_v3v3(edge2, nv3, nv2);
1725                                         sub_v3_v3v3(dv1, opco, nv2); /* abuse dv1 to have vertex in question at *origin* of triangle */
1726
1727                                         cross_v3_v3v3(d_nvect, edge2, edge1);
1728                                         /* n_mag = */ /* UNUSED */ normalize_v3(d_nvect);
1729                                         facedist = dot_v3v3(dv1, d_nvect);
1730                                         // so rules are
1731                                         //
1732
1733                                         if ((facedist > innerfacethickness) && (facedist < outerfacethickness)) {
1734                                                 if (isect_point_tri_prism_v3(opco, nv1, nv2, nv3) ) {
1735                                                         force_mag_norm =(float)exp(-ee*facedist);
1736                                                         if (facedist > outerfacethickness*ff)
1737                                                                 force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness);
1738                                                         *damp=ob->pd->pdef_sbdamp;
1739                                                         if (facedist > 0.0f) {
1740                                                                 *damp *= (1.0f - facedist/outerfacethickness);
1741                                                                 madd_v3_v3fl(outerforceaccu, d_nvect, force_mag_norm);
1742                                                                 deflected = 3;
1743
1744                                                         }
1745                                                         else {
1746                                                                 madd_v3_v3fl(innerforceaccu, d_nvect, force_mag_norm);
1747                                                                 if (deflected < 2) deflected = 2;
1748                                                         }
1749                                                         if ((mprevvert) && (*damp > 0.0f)) {
1750                                                                 choose_winner(ve, opco, nv1, nv2, nv3, vv1, vv2, vv3);
1751                                                                 add_v3_v3(avel, ve);
1752                                                                 cavel ++;
1753                                                         }
1754                                                         *intrusion += facedist;
1755                                                         ci++;
1756                                                 }
1757                                         }
1758
1759                                         mima++;
1760                                         vt++;
1761                                 }/* while a */
1762                         } /* if (ob->pd && ob->pd->deflect) */
1763                         BLI_ghashIterator_step(ihash);
1764                 }
1765         } /* while () */
1766
1767         if (deflected == 1) { // no face but 'outer' edge cylinder sees vert
1768                 force_mag_norm =(float)exp(-ee*mindistedge);
1769                 if (mindistedge > outerfacethickness*ff)
1770                         force_mag_norm =(float)force_mag_norm*fa*(mindistedge - outerfacethickness)*(mindistedge - outerfacethickness);
1771                 madd_v3_v3fl(force, coledge, force_mag_norm);
1772                 *damp=ob->pd->pdef_sbdamp;
1773                 if (mindistedge > 0.0f) {
1774                         *damp *= (1.0f - mindistedge/outerfacethickness);
1775                 }
1776
1777         }
1778         if (deflected == 2) { //  face inner detected
1779                 add_v3_v3(force, innerforceaccu);
1780         }
1781         if (deflected == 3) { //  face outer detected
1782                 add_v3_v3(force, outerforceaccu);
1783         }
1784
1785         BLI_ghashIterator_free(ihash);
1786         if (cavel) mul_v3_fl(avel, 1.0f/(float)cavel);
1787         copy_v3_v3(vel, avel);
1788         if (ci) *intrusion /= ci;
1789         if (deflected) {
1790                 normalize_v3_v3(facenormal, force);
1791         }
1792         return deflected;
1793 }
1794
1795
1796 /* sandbox to plug in various deflection algos */
1797 static int sb_deflect_face(Object *ob, float *actpos, float *facenormal, float *force, float *cf, float time, float *vel, float *intrusion)
1798 {
1799         float s_actpos[3];
1800         int deflected;
1801         copy_v3_v3(s_actpos, actpos);
1802         deflected= sb_detect_vertex_collisionCached(s_actpos, facenormal, cf, force, ob, time, vel, intrusion);
1803         //deflected= sb_detect_vertex_collisionCachedEx(s_actpos, facenormal, cf, force, ob, time, vel, intrusion);
1804         return(deflected);
1805 }
1806
1807 static void sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime))
1808 {
1809         SoftBody *sb= ob->soft; /* is supposed to be there */
1810         BodyPoint  *bp1, *bp2;
1811
1812         float dir[3], dvel[3];
1813         float distance, forcefactor, kd, absvel, projvel, kw;
1814
1815         /* prepare depending on which side of the spring we are on */
1816         if (bpi == bs->v1) {
1817                 bp1 = &sb->bpoint[bs->v1];
1818                 bp2 = &sb->bpoint[bs->v2];
1819         }
1820         else if (bpi == bs->v2) {
1821                 bp1 = &sb->bpoint[bs->v2];
1822                 bp2 = &sb->bpoint[bs->v1];
1823         }
1824         else {
1825                 /* TODO make this debug option */
1826                 /**/
1827                 printf("bodypoint <bpi> is not attached to spring  <*bs> --> sb_spring_force()\n");
1828                 return;
1829         }
1830
1831         /* do bp1 <--> bp2 elastic */
1832         sub_v3_v3v3(dir, bp1->pos, bp2->pos);
1833         distance = normalize_v3(dir);
1834         if (bs->len < distance)
1835                 iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
1836         else
1837                 iks  = 1.0f/(1.0f-sb->inpush)-1.0f ;/* inner spring constants function */
1838
1839         if (bs->len > 0.0f) /* check for degenerated springs */
1840                 forcefactor = iks/bs->len;
1841         else
1842                 forcefactor = iks;
1843         kw = (bp1->springweight+bp2->springweight)/2.0f;
1844         kw = kw * kw;
1845         kw = kw * kw;
1846         switch (bs->springtype) {
1847                 case SB_EDGE:
1848                 case SB_HANDLE:
1849                         forcefactor *=  kw;
1850                         break;
1851                 case SB_BEND:
1852                         forcefactor *=sb->secondspring*kw;
1853                         break;
1854                 case SB_STIFFQUAD:
1855                         forcefactor *=sb->shearstiff*sb->shearstiff* kw;
1856                         break;
1857                 default:
1858                         break;
1859         }
1860
1861
1862         madd_v3_v3fl(bp1->force, dir, (bs->len - distance) * forcefactor);
1863
1864         /* do bp1 <--> bp2 viscous */
1865         sub_v3_v3v3(dvel, bp1->vec, bp2->vec);
1866         kd = sb->infrict * sb_fric_force_scale(ob);
1867         absvel  = normalize_v3(dvel);
1868         projvel = dot_v3v3(dir, dvel);
1869         kd     *= absvel * projvel;
1870         madd_v3_v3fl(bp1->force, dir, -kd);
1871 }
1872
1873
1874 /* since this is definitely the most CPU consuming task here .. try to spread it */
1875 /* core function _softbody_calc_forces_slice_in_a_thread */
1876 /* result is int to be able to flag user break */
1877 static int _softbody_calc_forces_slice_in_a_thread(Scene *scene, Object *ob, float forcetime, float timenow, int ifirst, int ilast, int *UNUSED(ptr_to_break_func(void)), ListBase *effectors, int do_deflector, float fieldfactor, float windfactor)
1878 {
1879         float iks;
1880         int bb, do_selfcollision, do_springcollision, do_aero;
1881         int number_of_points_here = ilast - ifirst;
1882         SoftBody *sb= ob->soft; /* is supposed to be there */
1883         BodyPoint  *bp;
1884
1885         /* initialize */
1886         if (sb) {
1887                 /* check conditions for various options */
1888                 /* +++ could be done on object level to squeeze out the last bits of it */
1889                 do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
1890                 do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
1891                 do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
1892                 /* --- could be done on object level to squeeze out the last bits of it */
1893         }
1894         else {
1895                 printf("Error expected a SB here\n");
1896                 return (999);
1897         }
1898
1899 /* debugerin */
1900         if (sb->totpoint < ifirst) {
1901                 printf("Aye 998");
1902                 return (998);
1903         }
1904 /* debugerin */
1905
1906
1907         bp = &sb->bpoint[ifirst];
1908         for (bb=number_of_points_here; bb>0; bb--, bp++) {
1909                 /* clear forces  accumulator */
1910                 bp->force[0] = bp->force[1] = bp->force[2] = 0.0;
1911                 /* naive ball self collision */
1912                 /* needs to be done if goal snaps or not */
1913                 if (do_selfcollision) {
1914                         int attached;
1915                         BodyPoint   *obp;
1916                         BodySpring *bs;
1917                         int c, b;
1918                         float velcenter[3], dvel[3], def[3];
1919                         float distance;
1920                         float compare;
1921                         float bstune = sb->ballstiff;
1922
1923                         /* running in a slice we must not assume anything done with obp  neither alter the data of obp */
1924                         for (c=sb->totpoint, obp= sb->bpoint; c>0; c--, obp++) {
1925                                 compare = (obp->colball + bp->colball);
1926                                 sub_v3_v3v3(def, bp->pos, obp->pos);
1927                                 /* rather check the AABBoxes before ever calculating the real distance */
1928                                 /* mathematically it is completely nuts, but performance is pretty much (3) times faster */
1929                                 if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
1930                                 distance = normalize_v3(def);
1931                                 if (distance < compare ) {
1932                                         /* exclude body points attached with a spring */
1933                                         attached = 0;
1934                                         for (b=obp->nofsprings;b>0;b--) {
1935                                                 bs = sb->bspring + obp->springs[b-1];
1936                                                 if (( ilast-bb == bs->v2) || ( ilast-bb == bs->v1)) {
1937                                                         attached=1;
1938                                                         continue;}
1939                                         }
1940                                         if (!attached) {
1941                                                 float f = bstune / (distance) + bstune / (compare * compare) * distance - 2.0f * bstune / compare;
1942
1943                                                 mid_v3_v3v3(velcenter, bp->vec, obp->vec);
1944                                                 sub_v3_v3v3(dvel, velcenter, bp->vec);
1945                                                 mul_v3_fl(dvel, _final_mass(ob, bp));
1946
1947                                                 madd_v3_v3fl(bp->force, def, f * (1.0f - sb->balldamp));
1948                                                 madd_v3_v3fl(bp->force, dvel, sb->balldamp);
1949                                         }
1950                                 }
1951                         }
1952                 }
1953                 /* naive ball self collision done */
1954
1955                 if (_final_goal(ob, bp) < SOFTGOALSNAP) {  /* omit this bp when it snaps */
1956                         float auxvect[3];
1957                         float velgoal[3];
1958
1959                         /* do goal stuff */
1960                         if (ob->softflag & OB_SB_GOAL) {
1961                                 /* true elastic goal */
1962                                 float ks, kd;
1963                                 sub_v3_v3v3(auxvect, bp->pos, bp->origT);
1964                                 ks  = 1.0f / (1.0f - _final_goal(ob, bp) * sb->goalspring) - 1.0f;
1965                                 bp->force[0]+= -ks*(auxvect[0]);
1966                                 bp->force[1]+= -ks*(auxvect[1]);
1967                                 bp->force[2]+= -ks*(auxvect[2]);
1968
1969                                 /* calculate damping forces generated by goals*/
1970                                 sub_v3_v3v3(velgoal, bp->origS, bp->origE);
1971                                 kd =  sb->goalfrict * sb_fric_force_scale(ob);
1972                                 add_v3_v3v3(auxvect, velgoal, bp->vec);
1973
1974                                 if (forcetime > 0.0f) { /* make sure friction does not become rocket motor on time reversal */
1975                                         bp->force[0]-= kd * (auxvect[0]);
1976                                         bp->force[1]-= kd * (auxvect[1]);
1977                                         bp->force[2]-= kd * (auxvect[2]);
1978                                 }
1979                                 else {
1980                                         bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
1981                                         bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
1982                                         bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
1983                                 }
1984                         }
1985                         /* done goal stuff */
1986
1987                         /* gravitation */
1988                         if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
1989                                 float gravity[3];
1990                                 copy_v3_v3(gravity, scene->physics_settings.gravity);
1991                                 mul_v3_fl(gravity, sb_grav_force_scale(ob)*_final_mass(ob, bp)*sb->effector_weights->global_gravity); /* individual mass of node here */
1992                                 add_v3_v3(bp->force, gravity);
1993                         }
1994
1995                         /* particle field & vortex */
1996                         if (effectors) {
1997                                 EffectedPoint epoint;
1998                                 float kd;
1999                                 float force[3] = {0.0f, 0.0f, 0.0f};
2000                                 float speed[3] = {0.0f, 0.0f, 0.0f};
2001                                 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */
2002                                 pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint);
2003                                 BKE_effectors_apply(effectors, NULL, sb->effector_weights, &epoint, force, speed);
2004
2005                                 /* apply forcefield*/
2006                                 mul_v3_fl(force, fieldfactor* eval_sb_fric_force_scale);
2007                                 add_v3_v3(bp->force, force);
2008
2009                                 /* BP friction in moving media */
2010                                 kd= sb->mediafrict* eval_sb_fric_force_scale;
2011                                 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
2012                                 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
2013                                 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
2014                                 /* now we'll have nice centrifugal effect for vortex */
2015
2016                         }
2017                         else {
2018                                 /* BP friction in media (not) moving*/
2019                                 float kd = sb->mediafrict* sb_fric_force_scale(ob);
2020                                 /* assume it to be proportional to actual velocity */
2021                                 bp->force[0]-= bp->vec[0]*kd;
2022                                 bp->force[1]-= bp->vec[1]*kd;
2023                                 bp->force[2]-= bp->vec[2]*kd;
2024                                 /* friction in media done */
2025                         }
2026                         /* +++cached collision targets */
2027                         bp->choke = 0.0f;
2028                         bp->choke2 = 0.0f;
2029                         bp->loc_flag &= ~SBF_DOFUZZY;
2030                         if (do_deflector && !(bp->loc_flag & SBF_OUTOFCOLLISION) ) {
2031                                 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;
2032                                 float kd = 1.0f;
2033
2034                                 if (sb_deflect_face(ob, bp->pos, facenormal, defforce, &cf, timenow, vel, &intrusion)) {
2035                                         if (intrusion < 0.0f) {
2036                                                 sb->scratch->flag |= SBF_DOFUZZY;
2037                                                 bp->loc_flag |= SBF_DOFUZZY;
2038                                                 bp->choke = sb->choke*0.01f;
2039                                         }
2040
2041                                         sub_v3_v3v3(cfforce, bp->vec, vel);
2042                                         madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
2043
2044                                         madd_v3_v3fl(bp->force, defforce, kd);
2045                                 }
2046
2047                         }
2048                         /* ---cached collision targets */
2049
2050                         /* +++springs */
2051                         iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
2052                         if (ob->softflag & OB_SB_EDGES) {
2053                                 if (sb->bspring) { /* spring list exists at all ? */
2054                                         int b;
2055                                         BodySpring *bs;
2056                                         for (b=bp->nofsprings;b>0;b--) {
2057                                                 bs = sb->bspring + bp->springs[b-1];
2058                                                 if (do_springcollision || do_aero) {
2059                                                         add_v3_v3(bp->force, bs->ext_force);
2060                                                         if (bs->flag & BSF_INTERSECT)
2061                                                                 bp->choke = bs->cf;
2062
2063                                                 }
2064                                                 // sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float forcetime)
2065                                                 sb_spring_force(ob, ilast-bb, bs, iks, forcetime);
2066                                         }/* loop springs */
2067                                 }/* existing spring list */
2068                         }/*any edges*/
2069                         /* ---springs */
2070                 }/*omit on snap */
2071         }/*loop all bp's*/
2072         return 0; /*done fine*/
2073 }
2074
2075 static void *exec_softbody_calc_forces(void *data)
2076 {
2077         SB_thread_context *pctx = (SB_thread_context*)data;
2078         _softbody_calc_forces_slice_in_a_thread(pctx->scene, pctx->ob, pctx->forcetime, pctx->timenow, pctx->ifirst, pctx->ilast, NULL, pctx->effectors, pctx->do_deflector, pctx->fieldfactor, pctx->windfactor);
2079         return NULL;
2080 }
2081
2082 static void sb_cf_threads_run(Scene *scene, Object *ob, float forcetime, float timenow, int totpoint, int *UNUSED(ptr_to_break_func(void)), struct ListBase *effectors, int do_deflector, float fieldfactor, float windfactor)
2083 {
2084         ListBase threads;
2085         SB_thread_context *sb_threads;
2086         int i, totthread, left, dec;
2087         int lowpoints =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */
2088
2089         /* figure the number of threads while preventing pretty pointless threading overhead */
2090         totthread= BKE_scene_num_threads(scene);
2091         /* what if we got zillions of CPUs running but less to spread*/
2092         while ((totpoint/totthread < lowpoints) && (totthread > 1)) {
2093                 totthread--;
2094         }
2095
2096         /* printf("sb_cf_threads_run spawning %d threads\n", totthread); */
2097
2098         sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBThread");
2099         memset(sb_threads, 0, sizeof(SB_thread_context)*totthread);
2100         left = totpoint;
2101         dec = totpoint/totthread +1;
2102         for (i=0; i<totthread; i++) {
2103                 sb_threads[i].scene = scene;
2104                 sb_threads[i].ob = ob;
2105                 sb_threads[i].forcetime = forcetime;
2106                 sb_threads[i].timenow = timenow;
2107                 sb_threads[i].ilast   = left;
2108                 left = left - dec;
2109                 if (left >0) {
2110                         sb_threads[i].ifirst  = left;
2111                 }
2112                 else
2113                         sb_threads[i].ifirst  = 0;
2114                 sb_threads[i].effectors = effectors;
2115                 sb_threads[i].do_deflector = do_deflector;
2116                 sb_threads[i].fieldfactor = fieldfactor;
2117                 sb_threads[i].windfactor  = windfactor;
2118                 sb_threads[i].nr= i;
2119                 sb_threads[i].tot= totthread;
2120         }
2121
2122
2123         if (totthread > 1) {
2124                 BLI_threadpool_init(&threads, exec_softbody_calc_forces, totthread);
2125
2126                 for (i=0; i<totthread; i++)
2127                         BLI_threadpool_insert(&threads, &sb_threads[i]);
2128
2129                 BLI_threadpool_end(&threads);
2130         }
2131         else
2132                 exec_softbody_calc_forces(&sb_threads[0]);
2133         /* clean up */
2134         MEM_freeN(sb_threads);
2135 }
2136
2137 static void softbody_calc_forcesEx(struct Depsgraph *depsgraph, Scene *scene, Object *ob, float forcetime, float timenow)
2138 {
2139         /* rule we never alter free variables :bp->vec bp->pos in here !
2140          * this will ruin adaptive stepsize AKA heun! (BM)
2141          */
2142         SoftBody *sb= ob->soft; /* is supposed to be there */
2143         /*BodyPoint *bproot;*/ /* UNUSED */
2144         /* float gravity; */ /* UNUSED */
2145         /* float iks; */
2146         float fieldfactor = -1.0f, windfactor  = 0.25;
2147         int   do_deflector /*, do_selfcollision*/, do_springcollision, do_aero;
2148
2149         /* gravity = sb->grav * sb_grav_force_scale(ob); */ /* UNUSED */
2150
2151         /* check conditions for various options */
2152         do_deflector= query_external_colliders(depsgraph, sb->collision_group);
2153         /* do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); */ /* UNUSED */
2154         do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
2155         do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
2156
2157         /* iks  = 1.0f/(1.0f-sb->inspring)-1.0f; */ /* inner spring constants function */ /* UNUSED */
2158         /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */
2159
2160         if (do_springcollision || do_aero)
2161                 sb_sfesf_threads_run(depsgraph, scene, ob, timenow, sb->totspring, NULL);
2162
2163         /* after spring scan because it uses Effoctors too */
2164         ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, sb->effector_weights);
2165
2166         if (do_deflector) {
2167                 float defforce[3];
2168                 do_deflector = sb_detect_aabb_collisionCached(defforce, ob, timenow);
2169         }
2170
2171         sb_cf_threads_run(scene, ob, forcetime, timenow, sb->totpoint, NULL, effectors, do_deflector, fieldfactor, windfactor);
2172
2173         /* finally add forces caused by face collision */
2174         if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob, timenow);
2175
2176         /* finish matrix and solve */
2177         BKE_effectors_free(effectors);
2178 }
2179
2180
2181 static void softbody_calc_forces(struct Depsgraph *depsgraph, Scene *scene, Object *ob, float forcetime, float timenow)
2182 {
2183         /* redirection to the new threaded Version */
2184         if (!(G.debug_value & 0x10)) { // 16
2185                 softbody_calc_forcesEx(depsgraph, scene, ob, forcetime, timenow);
2186                 return;
2187         }
2188         else {
2189                 /* so the following will die  */
2190                 /* |||||||||||||||||||||||||| */
2191                 /* VVVVVVVVVVVVVVVVVVVVVVVVVV */
2192                 /*backward compatibility note:
2193                 fixing bug [17428] which forces adaptive step size to tiny steps
2194                 in some situations
2195                 .. keeping G.debug_value==17 0x11 option for old files 'needing' the bug*/
2196
2197                 /* rule we never alter free variables :bp->vec bp->pos in here !
2198                  * this will ruin adaptive stepsize AKA heun! (BM)
2199                  */
2200                 SoftBody *sb= ob->soft; /* is supposed to be there */
2201                 BodyPoint  *bp;
2202                 /* BodyPoint *bproot; */ /* UNUSED */
2203                 BodySpring *bs;
2204                 float iks, ks, kd, gravity[3] = {0.0f, 0.0f, 0.0f};
2205                 float fieldfactor = -1.0f, windfactor  = 0.25f;
2206                 float tune = sb->ballstiff;
2207                 int do_deflector, do_selfcollision, do_springcollision, do_aero;
2208
2209                 if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
2210                         copy_v3_v3(gravity, scene->physics_settings.gravity);
2211                         mul_v3_fl(gravity, sb_grav_force_scale(ob)*sb->effector_weights->global_gravity);
2212                 }
2213
2214                 /* check conditions for various options */
2215                 do_deflector= query_external_colliders(depsgraph, sb->collision_group);
2216                 do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
2217                 do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL);
2218                 do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES));
2219
2220                 iks  = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
2221                 /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */
2222
2223                 if (do_springcollision || do_aero)  scan_for_ext_spring_forces(depsgraph, scene, ob, timenow);
2224                 /* after spring scan because it uses Effoctors too */
2225                 ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, ob->soft->effector_weights);
2226
2227                 if (do_deflector) {
2228                         float defforce[3];
2229                         do_deflector = sb_detect_aabb_collisionCached(defforce, ob, timenow);
2230                 }
2231
2232                 bp = sb->bpoint;
2233                 for (int a = sb->totpoint; a > 0; a--, bp++) {
2234                         /* clear forces  accumulator */
2235                         bp->force[0] = bp->force[1] = bp->force[2] = 0.0;
2236
2237                         /* naive ball self collision */
2238                         /* needs to be done if goal snaps or not */
2239                         if (do_selfcollision) {
2240                                 int attached;
2241                                 BodyPoint   *obp;
2242                                 int c, b;
2243                                 float velcenter[3], dvel[3], def[3];
2244                                 float distance;
2245                                 float compare;
2246
2247                                 for (c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) {
2248
2249                                         //if ((bp->octantflag & obp->octantflag) == 0) continue;
2250
2251                                         compare = (obp->colball + bp->colball);
2252                                         sub_v3_v3v3(def, bp->pos, obp->pos);
2253
2254                                         /* rather check the AABBoxes before ever calculating the real distance */
2255                                         /* mathematically it is completely nuts, but performance is pretty much (3) times faster */
2256                                         if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
2257
2258                                         distance = normalize_v3(def);
2259                                         if (distance < compare ) {
2260                                                 /* exclude body points attached with a spring */
2261                                                 attached = 0;
2262                                                 for (b=obp->nofsprings;b>0;b--) {
2263                                                         bs = sb->bspring + obp->springs[b-1];
2264                                                         if (( sb->totpoint-a == bs->v2) || ( sb->totpoint-a == bs->v1)) {
2265                                                                 attached=1;
2266                                                                 continue;}
2267                                                 }
2268                                                 if (!attached) {
2269                                                         float f = tune / (distance) + tune / (compare * compare) * distance - 2.0f * tune/compare;
2270
2271                                                         mid_v3_v3v3(velcenter, bp->vec, obp->vec);
2272                                                         sub_v3_v3v3(dvel, velcenter, bp->vec);
2273                                                         mul_v3_fl(dvel, _final_mass(ob, bp));
2274
2275                                                         madd_v3_v3fl(bp->force, def, f * (1.0f - sb->balldamp));
2276                                                         madd_v3_v3fl(bp->force, dvel, sb->balldamp);
2277
2278                                                         /* exploit force(a, b) == -force(b, a) part2/2 */
2279                                                         sub_v3_v3v3(dvel, velcenter, obp->vec);
2280                                                         mul_v3_fl(dvel, (_final_mass(ob, bp)+_final_mass(ob, obp))/2.0f);
2281
2282                                                         madd_v3_v3fl(obp->force, dvel, sb->balldamp);
2283                                                         madd_v3_v3fl(obp->force, def, -f * (1.0f - sb->balldamp));
2284                                                 }
2285                                         }
2286                                 }
2287                         }
2288                         /* naive ball self collision done */
2289
2290                         if (_final_goal(ob, bp) < SOFTGOALSNAP) {  /* omit this bp when it snaps */
2291                                 float auxvect[3];
2292                                 float velgoal[3];
2293
2294                                 /* do goal stuff */
2295                                 if (ob->softflag & OB_SB_GOAL) {
2296                                         /* true elastic goal */
2297                                         sub_v3_v3v3(auxvect, bp->pos, bp->origT);
2298                                         ks  = 1.0f / (1.0f- _final_goal(ob, bp) * sb->goalspring) - 1.0f;
2299                                         bp->force[0]+= -ks*(auxvect[0]);
2300                                         bp->force[1]+= -ks*(auxvect[1]);
2301                                         bp->force[2]+= -ks*(auxvect[2]);
2302
2303                                         /* calculate damping forces generated by goals*/
2304                                         sub_v3_v3v3(velgoal, bp->origS, bp->origE);
2305                                         kd = sb->goalfrict * sb_fric_force_scale(ob);
2306                                         add_v3_v3v3(auxvect, velgoal, bp->vec);
2307
2308                                         if (forcetime > 0.0f) { /* make sure friction does not become rocket motor on time reversal */
2309                                                 bp->force[0]-= kd * (auxvect[0]);
2310                                                 bp->force[1]-= kd * (auxvect[1]);
2311                                                 bp->force[2]-= kd * (auxvect[2]);
2312
2313                                         }
2314                                         else {
2315                                                 bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
2316                                                 bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
2317                                                 bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
2318                                         }
2319                                 }
2320                                 /* done goal stuff */
2321
2322
2323                                 /* gravitation */
2324                                 madd_v3_v3fl(bp->force, gravity, _final_mass(ob, bp)); /* individual mass of node here */
2325
2326
2327                                 /* particle field & vortex */
2328                                 if (effectors) {
2329                                         EffectedPoint epoint;
2330                                         float force[3] = {0.0f, 0.0f, 0.0f};
2331                                         float speed[3] = {0.0f, 0.0f, 0.0f};
2332                                         float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */
2333                                         pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint);
2334                                         BKE_effectors_apply(effectors, NULL, sb->effector_weights, &epoint, force, speed);
2335
2336                                         /* apply forcefield*/
2337                                         mul_v3_fl(force, fieldfactor* eval_sb_fric_force_scale);
2338                                         add_v3_v3(bp->force, force);
2339
2340                                         /* BP friction in moving media */
2341                                         kd= sb->mediafrict* eval_sb_fric_force_scale;
2342                                         bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
2343                                         bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
2344                                         bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
2345                                         /* now we'll have nice centrifugal effect for vortex */
2346
2347                                 }
2348                                 else {
2349                                         /* BP friction in media (not) moving*/
2350                                         kd= sb->mediafrict* sb_fric_force_scale(ob);
2351                                         /* assume it to be proportional to actual velocity */
2352                                         bp->force[0]-= bp->vec[0]*kd;
2353                                         bp->force[1]-= bp->vec[1]*kd;
2354                                         bp->force[2]-= bp->vec[2]*kd;
2355                                         /* friction in media done */
2356                                 }
2357                                 /* +++cached collision targets */
2358                                 bp->choke = 0.0f;
2359                                 bp->choke2 = 0.0f;
2360                                 bp->loc_flag &= ~SBF_DOFUZZY;
2361                                 if (do_deflector) {
2362                                         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;
2363                                         kd = 1.0f;
2364
2365                                         if (sb_deflect_face(ob, bp->pos, facenormal, defforce, &cf, timenow, vel, &intrusion)) {
2366                                                 if (intrusion < 0.0f) {
2367                                                         if (G.debug_value & 0x01) { // 17 we did check for bit 0x10 before
2368                                                                 /* fixing bug [17428] this forces adaptive step size to tiny steps
2369                                                                  * in some situations .. keeping G.debug_value==17 option for old files 'needing' the bug
2370                                                                  */
2371                                                                 /* bjornmose:  uugh.. what an evil hack
2372                                                                  * violation of the 'don't touch bp->pos in here' rule
2373                                                                  * but works nice, like this-->
2374                                                                  * we predict the solution being out of the collider
2375                                                                  * in heun step No1 and leave the heun step No2 adapt to it
2376                                                                  * so we kind of introduced a implicit solver for this case
2377                                                                  */
2378                                                                 madd_v3_v3fl(bp->pos, facenormal, -intrusion);
2379                                                         }
2380                                                         else {
2381
2382                                                                 sub_v3_v3v3(cfforce, bp->vec, vel);
2383                                                                 madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
2384                                                         }
2385
2386
2387                                                         sb->scratch->flag |= SBF_DOFUZZY;
2388                                                         bp->loc_flag |= SBF_DOFUZZY;
2389                                                         bp->choke = sb->choke*0.01f;
2390                                                 }
2391                                                 else {
2392                                                         sub_v3_v3v3(cfforce, bp->vec, vel);
2393                                                         madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
2394                                                 }
2395                                                 madd_v3_v3fl(bp->force, defforce, kd);
2396
2397                                         }
2398
2399                                 }
2400                                 /* ---cached collision targets */
2401
2402                                 /* +++springs */
2403                                 if (ob->softflag & OB_SB_EDGES) {
2404                                         if (sb->bspring) { /* spring list exists at all ? */
2405                                                 for (int b = bp->nofsprings; b > 0; b--) {
2406                                                         bs = sb->bspring + bp->springs[b-1];
2407                                                         if (do_springcollision || do_aero) {
2408                                                                 add_v3_v3(bp->force, bs->ext_force);
2409                                                                 if (bs->flag & BSF_INTERSECT)
2410                                                                         bp->choke = bs->cf;
2411
2412                                                         }
2413                                                         // sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float forcetime)
2414                                                         // rather remove nl_falgs from code .. will make things a lot cleaner
2415                                                         sb_spring_force(ob, sb->totpoint-a, bs, iks, forcetime);
2416                                                 }/* loop springs */
2417                                         }/* existing spring list */
2418                                 }/*any edges*/
2419                                 /* ---springs */
2420                         }/*omit on snap */
2421                 }/*loop all bp's*/
2422
2423
2424                 /* finally add forces caused by face collision */
2425                 if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob, timenow);
2426                 BKE_effectors_free(effectors);
2427         }
2428 }
2429
2430 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags)
2431 {
2432         /* time evolution */
2433         /* actually does an explicit euler step mode == 0 */
2434         /* or heun ~ 2nd order runge-kutta steps, mode 1, 2 */
2435         SoftBody *sb= ob->soft; /* is supposed to be there */
2436         BodyPoint *bp;
2437         float dx[3] = {0}, dv[3], aabbmin[3], aabbmax[3], cm[3] = {0.0f, 0.0f, 0.0f};
2438         float timeovermass/*, freezeloc=0.00001f, freezeforce=0.00000000001f*/;
2439         float maxerrpos= 0.0f, maxerrvel = 0.0f;
2440         int a, fuzzy=0;
2441
2442         forcetime *= sb_time_scale(ob);
2443
2444         aabbmin[0]=aabbmin[1]=aabbmin[2] = 1e20f;
2445         aabbmax[0]=aabbmax[1]=aabbmax[2] = -1e20f;
2446
2447         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
2448 /* now we have individual masses   */
2449 /* claim a minimum mass for vertex */
2450                 if (_final_mass(ob, bp) > 0.009999f) timeovermass = forcetime/_final_mass(ob, bp);
2451                 else timeovermass = forcetime/0.009999f;
2452
2453
2454                 if (_final_goal(ob, bp) < SOFTGOALSNAP) {
2455                         /* this makes t~ = t */
2456                         if (mid_flags & MID_PRESERVE) copy_v3_v3(dx, bp->vec);
2457
2458                         /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces  + more forces*/
2459                         /* the ( ... )' operator denotes derivate respective time */
2460                         /* the euler step for velocity then becomes */
2461                         /* v(t + dt) = v(t) + a(t) * dt */
2462                         mul_v3_fl(bp->force, timeovermass);/* individual mass of node here */
2463                         /* some nasty if's to have heun in here too */
2464                         copy_v3_v3(dv, bp->force);
2465
2466                         if (mode == 1) {
2467                                 copy_v3_v3(bp->prevvec, bp->vec);
2468                                 copy_v3_v3(bp->prevdv, dv);
2469                         }
2470
2471                         if (mode ==2) {
2472                                 /* be optimistic and execute step */
2473                                 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]);
2474                                 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]);
2475                                 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]);
2476                                 /* compare euler to heun to estimate error for step sizing */
2477                                 maxerrvel = max_ff(maxerrvel, fabsf(dv[0] - bp->prevdv[0]));
2478                                 maxerrvel = max_ff(maxerrvel, fabsf(dv[1] - bp->prevdv[1]));
2479                                 maxerrvel = max_ff(maxerrvel, fabsf(dv[2] - bp->prevdv[2]));
2480                         }
2481                         else { add_v3_v3(bp->vec, bp->force); }
2482
2483                         /* this makes t~ = t+dt */
2484                         if (!(mid_flags & MID_PRESERVE)) copy_v3_v3(dx, bp->vec);
2485
2486                         /* so here is (x)'= v(elocity) */
2487                         /* the euler step for location then becomes */
2488                         /* x(t + dt) = x(t) + v(t~) * dt */
2489                         mul_v3_fl(dx, forcetime);
2490
2491                         /* again some nasty if's to have heun in here too */
2492                         if (mode ==1) {
2493                                 copy_v3_v3(bp->prevpos, bp->pos);
2494                                 copy_v3_v3(bp->prevdx, dx);
2495                         }
2496
2497                         if (mode ==2) {
2498                                 bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]);
2499                                 bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]);
2500                                 bp->pos[2] = bp->prevpos[2] + 0.5f * ( dx[2] + bp->prevdx[2]);
2501                                 maxerrpos = max_ff(maxerrpos, fabsf(dx[0] - bp->prevdx[0]));
2502                                 maxerrpos = max_ff(maxerrpos, fabsf(dx[1] - bp->prevdx[1]));
2503                                 maxerrpos = max_ff(maxerrpos, fabsf(dx[2] - bp->prevdx[2]));
2504
2505                                 /* bp->choke is set when we need to pull a vertex or edge out of the collider.
2506                                  * the collider object signals to get out by pushing hard. on the other hand
2507                                  * we don't want to end up in deep space so we add some <viscosity>
2508                                  * to balance that out */
2509                                 if (bp->choke2 > 0.0f) {
2510                                         mul_v3_fl(bp->vec, (1.0f - bp->choke2));
2511                                 }
2512                                 if (bp->choke > 0.0f) {
2513                                         mul_v3_fl(bp->vec, (1.0f - bp->choke));
2514                                 }
2515
2516                         }
2517                         else { add_v3_v3(bp->pos, dx);}
2518                 }/*snap*/
2519                 /* so while we are looping BPs anyway do statistics on the fly */
2520                 minmax_v3v3_v3(aabbmin, aabbmax, bp->pos);
2521                 if (bp->loc_flag & SBF_DOFUZZY) fuzzy =1;
2522         } /*for*/
2523
2524         if (sb->totpoint) mul_v3_fl(cm, 1.0f/sb->totpoint);
2525         if (sb->scratch) {
2526                 copy_v3_v3(sb->scratch->aabbmin, aabbmin);
2527                 copy_v3_v3(sb->scratch->aabbmax, aabbmax);
2528         }
2529
2530         if (err) { /* so step size will be controlled by biggest difference in slope */
2531                 if (sb->solverflags & SBSO_OLDERR)
2532                         *err = max_ff(maxerrpos, maxerrvel);
2533                 else
2534                         *err = maxerrpos;
2535                 //printf("EP %f EV %f\n", maxerrpos, maxerrvel);
2536                 if (fuzzy) {
2537                         *err /= sb->fuzzyness;
2538                 }
2539         }
2540 }
2541
2542 /* used by heun when it overshoots */
2543 static void softbody_restore_prev_step(Object *ob)
2544 {
2545         SoftBody *sb= ob->soft; /* is supposed to be there*/
2546         BodyPoint *bp;
2547         int a;
2548
2549         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
2550                 copy_v3_v3(bp->vec, bp->prevvec);
2551                 copy_v3_v3(bp->pos, bp->prevpos);
2552         }
2553 }
2554
2555 /* care for bodypoints taken out of the 'ordinary' solver step
2556  * because they are screwed to goal by bolts
2557  * they just need to move along with the goal in time
2558  * we need to adjust them on sub frame timing in solver
2559  * so now when frame is done .. put 'em to the position at the end of frame
2560  */
2561 static void softbody_apply_goalsnap(Object *ob)
2562 {
2563         SoftBody *sb= ob->soft; /* is supposed to be there */
2564         BodyPoint *bp;
2565         int a;
2566
2567         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
2568                 if (_final_goal(ob, bp) >= SOFTGOALSNAP) {
2569                         copy_v3_v3(bp->prevpos, bp->pos);
2570                         copy_v3_v3(bp->pos, bp->origT);
2571                 }
2572         }
2573 }
2574
2575
2576 static void apply_spring_memory(Object *ob)
2577 {
2578         SoftBody *sb = ob->soft;
2579         BodySpring *bs;
2580         BodyPoint *bp1, *bp2;
2581         int a;
2582         float b, l, r;
2583
2584         if (sb && sb->totspring) {
2585                 b = sb->plastic;
2586                 for (a=0; a<sb->totspring; a++) {
2587                         bs  = &sb->bspring[a];
2588                         bp1 =&sb->bpoint[bs->v1];
2589                         bp2 =&sb->bpoint[bs->v2];
2590                         l = len_v3v3(bp1->pos, bp2->pos);
2591                         r = bs->len/l;
2592                         if (( r > 1.05f) || (r < 0.95f)) {
2593                                 bs->len = ((100.0f - b) * bs->len  + b*l)/100.0f;
2594                         }
2595                 }
2596         }
2597 }
2598
2599 /* expects full initialized softbody */
2600 static void interpolate_exciter(Object *ob, int timescale, int time)
2601 {
2602         SoftBody *sb= ob->soft;
2603         BodyPoint *bp;
2604         float f;
2605         int a;
2606
2607         f = (float)time/(float)timescale;
2608
2609         for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
2610                 bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]);
2611                 bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]);
2612                 bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]);
2613                 if (_final_goal(ob, bp) >= SOFTGOALSNAP) {
2614                         bp->vec[0] = bp->origE[0] - bp->origS[0];
2615                         bp->vec[1] = bp->origE[1] - bp->origS[1];
2616                         bp->vec[2] = bp->origE[2] - bp->origS[2];
2617                 }
2618         }
2619
2620 }
2621
2622
2623 /* ************ convertors ********** */
2624
2625 /* for each object type we need;
2626  * - xxxx_to_softbody(Object *ob)      : a full (new) copy, creates SB geometry
2627  */
2628
2629 /* Resetting a Mesh SB object's springs */
2630 /* Spring length are caculted from'raw' mesh vertices that are NOT altered by modifier stack. */
2631 static void springs_from_mesh(Object *ob)
2632 {
2633         SoftBody *sb;
2634         Mesh *me= ob->data;
2635         BodyPoint *bp;
2636         int a;
2637         float scale =1.0f;
2638
2639         sb= ob->soft;
2640         if (me && sb) {
2641                 /* using bp->origS as a container for spring calculations here
2642                  * will be overwritten sbObjectStep() to receive
2643                  * actual modifier stack positions
2644                  */
2645                 if (me->totvert) {
2646                         bp= ob->soft->bpoint;
2647                         for (a=0; a<me->totvert; a++, bp++) {
2648                                 copy_v3_v3(bp->origS, me->mvert[a].co);
2649                                 mul_m4_v3(ob->obmat, bp->origS);
2650                         }
2651
2652                 }
2653                 /* recalculate spring length for meshes here */
2654                 /* public version shrink to fit */
2655                 if (sb->springpreload != 0 ) {
2656                         scale = sb->springpreload / 100.0f;
2657                 }
2658                 for (a=0; a<sb->totspring; a++) {
2659                         BodySpring *bs = &sb->bspring[a];
2660                         bs->len= scale*len_v3v3(sb->bpoint[bs->v1].origS, sb->bpoint[bs->v2].origS);
2661                 }
2662         }
2663 }
2664
2665
2666
2667
2668 /* makes totally fresh start situation */
2669 static void mesh_to_softbody(Scene *scene, Object *ob)
2670 {
2671         SoftBody *sb;
2672         Mesh *me= ob->data;
2673         MEdge *medge= me->medge;
2674         BodyPoint *bp;
2675         BodySpring *bs;
2676         int a, totedge;
2677         int defgroup_index, defgroup_index_mass, defgroup_index_spring;
2678
2679         if (ob->softflag & OB_SB_EDGES) totedge= me->totedge;
2680         else totedge= 0;
2681
2682         /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
2683         renew_softbody(scene, ob, me->totvert, totedge);
2684
2685         /* we always make body points */
2686         sb = ob->soft;
2687         bp= sb->bpoint;
2688
2689         defgroup_index        = me->dvert ? (sb->vertgroup - 1) : -1;
2690         defgroup_index_mass   = me->dvert ? defgroup_name_index(ob, sb->namedVG_Mass) : -1;
2691         defgroup_index_spring = me->dvert ? defgroup_name_index(ob, sb->namedVG_Spring_K) : -1;
2692
2693         for (a=0; a<me->totvert; a++, bp++) {
2694                 /* get scalar values needed  *per vertex* from vertex group functions,
2695                  * so we can *paint* them nicly ..
2696                  * they are normalized [0.0..1.0] so may be we need amplitude for scale
2697                  * which can be done by caller but still .. i'd like it to go this way
2698                  */
2699
2700                 if (ob->softflag & OB_SB_GOAL) {
2701                         BLI_assert(bp->goal == sb->defgoal);
2702                 }
2703                 if ((ob->softflag & OB_SB_GOAL) && (defgroup_index != -1)) {
2704                         bp->goal *= defvert_find_weight(&me->dvert[a], defgroup_index);
2705                 }
2706
2707                 /* to proof the concept
2708                  * this enables per vertex *mass painting*
2709                  */
2710
2711                 if (defgroup_index_mass != -1) {
2712                         bp->mass *= defvert_find_weight(&me->dvert[a], defgroup_index_mass);
2713                 }
2714
2715                 if (defgroup_index_spring != -1) {
2716                         bp->springweight *= defvert_find_weight(&me->dvert[a], defgroup_index_spring);
2717                 }
2718         }
2719
2720         /* but we only optionally add body edge springs */
2721         if (ob->softflag & OB_SB_EDGES) {
2722                 if (medge) {
2723                         bs= sb->bspring;
2724                         for (a=me->totedge; a>0; a--, medge++, bs++) {
2725                                 bs->v1= medge->v1;
2726                                 bs->v2= medge->v2;
2727                                 bs->springtype=SB_EDGE;
2728                         }
2729
2730
2731                         /* insert *diagonal* springs in quads if desired */
2732                         if (ob->softflag & OB_SB_QUADS) {
2733                                 add_mesh_quad_diag_springs(ob);
2734                         }
2735
2736                         build_bps_springlist(ob); /* scan for springs attached to bodypoints ONCE */
2737                         /* insert *other second order* springs if desired */
2738                         if (sb->secondspring > 0.0000001f) {
2739                                 add_2nd_order_springs(ob, sb->secondspring); /* exploits the first run of build_bps_springlist(ob);*/
2740                                 build_bps_springlist(ob); /* yes we need to do it again*/
2741                         }
2742                         springs_from_mesh(ob); /* write the 'rest'-length of the springs */
2743                         if (ob->softflag & OB_SB_SELF) {calculate_collision_balls(ob);}
2744
2745                 }
2746
2747         }
2748
2749 }
2750
2751 static void mesh_faces_to_scratch(Object *ob)
2752 {
2753         SoftBody *sb= ob->soft;
2754         const Mesh *me = ob->data;
2755         MLoopTri *looptri, *lt;
2756         BodyFace *bodyface;
2757         int a;
2758         /* alloc and copy faces*/
2759
2760         sb->scratch->totface = poly_to_tri_count(me->totpoly, me->totloop);
2761         looptri = lt = MEM_mallocN(sizeof(*looptri) * sb->scratch->totface, __func__);
2762         BKE_mesh_recalc_looptri(
2763                 me->mloop, me->mpoly,
2764                 me->mvert,
2765                 me->totloop, me->totpoly,
2766                 looptri);
2767
2768         bodyface = sb->scratch->bodyface = MEM_mallocN(sizeof(BodyFace) * sb->scratch->totface, "SB_body_Faces");
2769
2770         for (a = 0; a < sb->scratch->totface; a++, lt++, bodyface++) {
2771                 bodyface->v1 = me->mloop[lt->tri[0]].v;
2772                 bodyface->v2 = me->mloop[lt->tri[1]].v;
2773                 bodyface->v3 = me->mloop[lt->tri[2]].v;
2774                 zero_v3(bodyface->ext_force);
2775                 bodyface->ext_force[0] = bodyface->ext_force[1] = bodyface->ext_force[2] = 0.0f;
2776                 bodyface->flag = 0;
2777         }
2778
2779         MEM_freeN(looptri);
2780 }
2781 static void reference_to_scratch(Object *ob)
2782 {
2783         SoftBody *sb= ob->soft;
2784         ReferenceVert *rp;
2785         BodyPoint     *bp;
2786         float accu_pos[3] ={0.f, 0.f, 0.f};
2787         float accu_mass = 0.f;
2788         int a;
2789
2790         sb->scratch->Ref.ivert = MEM_mallocN(sizeof(ReferenceVert)*sb->totpoint, "SB_Reference");
2791         bp= ob->soft->bpoint;
2792         rp= sb->scratch->Ref.ivert;
2793         for (a=0; a<sb->totpoint; a++, rp++, bp++) {
2794                 copy_v3_v3(rp->pos, bp->pos);
2795                 add_v3_v3(accu_pos, bp->pos);
2796                 accu_mass += _final_mass(ob, bp);
2797         }
2798         mul_v3_fl(accu_pos, 1.0f/accu_mass);
2799         copy_v3_v3(sb->scratch->Ref.com, accu_pos);
2800         /* printf("reference_to_scratch\n"); */
2801 }
2802
2803 /*
2804  * helper function to get proper spring length
2805  * when object is rescaled
2806  */
2807 static float globallen(float *v1, float *v2, Object *ob)
2808 {
2809         float p1[3], p2[3];
2810         copy_v3_v3(p1, v1);
2811         mul_m4_v3(ob->obmat, p1);
2812         copy_v3_v3(p2, v2);
2813         mul_m4_v3(ob->obmat, p2);
2814         return len_v3v3(p1, p2);
2815 }
2816
2817 static void makelatticesprings(Lattice *lt,     BodySpring *bs, int dostiff, Object *ob)
2818 {
2819         BPoint *bp=lt->def, *bpu;
2820         int u, v, w, dv, dw, bpc=0, bpuc;
2821
2822         dv= lt->pntsu;
2823         dw= dv*lt->pntsv;
2824
2825         for (w=0; w<lt->pntsw; w++) {
2826
2827                 for (v=0; v<lt->pntsv; v++) {
2828
2829                         for (u=0, bpuc=0, bpu=NULL; u<lt->pntsu; u++, bp++, bpc++) {
2830
2831                                 if (w) {
2832                                         bs->v1 = bpc;
2833                                         bs->v2 = bpc-dw;
2834                                         bs->springtype=SB_EDGE;
2835