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