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