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