1 /*
2 * This program is free software; you can redistribute it and/or
3 * modify it under the terms of the GNU General Public License
4 * as published by the Free Software Foundation; either version 2
5 * of the License, or (at your option) any later version.
6 *
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
11 *
12 * You should have received a copy of the GNU General Public License
13 * along with this program; if not, write to the Free Software Foundation,
14 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15 *
16 * The Original Code is Copyright (C) Blender Foundation
17 * All rights reserved.
18 */
19
20 /** \file
21 * \ingroup bke
22 */
23
24 /**
25 * variables on the UI for now
26 * <pre>
27 * float mediafrict; friction to env
28 * float nodemass; softbody mass of *vertex*
29 * float grav; softbody amount of gravitation to apply
30 *
31 * float goalspring; softbody goal springs
32 * float goalfrict; softbody goal springs friction
33 * float mingoal; quick limits for goal
34 * float maxgoal;
35 *
36 * float inspring; softbody inner springs
37 * float infrict; softbody inner springs friction
38 * </pre>
39 */
40
41 #include <math.h>
42 #include <stdlib.h>
43 #include <string.h>
44
45 #include "CLG_log.h"
46
47 #include "MEM_guardedalloc.h"
48
49 /* types */
50 #include "DNA_collection_types.h"
51 #include "DNA_curve_types.h"
52 #include "DNA_lattice_types.h"
53 #include "DNA_mesh_types.h"
54 #include "DNA_meshdata_types.h"
55 #include "DNA_object_types.h"
56 #include "DNA_scene_types.h"
57
58 #include "BLI_ghash.h"
59 #include "BLI_listbase.h"
60 #include "BLI_math.h"
61 #include "BLI_threads.h"
62 #include "BLI_utildefines.h"
63
64 #include "BKE_collection.h"
65 #include "BKE_collision.h"
66 #include "BKE_curve.h"
67 #include "BKE_deform.h"
68 #include "BKE_effect.h"
69 #include "BKE_global.h"
70 #include "BKE_layer.h"
71 #include "BKE_mesh.h"
72 #include "BKE_modifier.h"
73 #include "BKE_pointcache.h"
74 #include "BKE_scene.h"
75 #include "BKE_softbody.h"
76
77 #include "DEG_depsgraph.h"
78 #include "DEG_depsgraph_query.h"
79
80 #include "PIL_time.h"
81
82 static CLG_LogRef LOG = {"bke.softbody"};
83
84 /* callbacks for errors and interrupts and some goo */
85 static int (*SB_localInterruptCallBack)(void) = NULL;
86
87 /* ********** soft body engine ******* */
88
89 typedef enum { SB_EDGE = 1, SB_BEND = 2, SB_STIFFQUAD = 3, SB_HANDLE = 4 } type_spring;
90
91 typedef struct BodySpring {
92 int v1, v2;
93 float len, cf, load;
94 float ext_force[3]; /* edges colliding and sailing */
95 type_spring springtype;
96 short flag;
97 } BodySpring;
98
99 typedef struct BodyFace {
100 int v1, v2, v3;
101 float ext_force[3]; /* faces colliding */
102 short flag;
103 } BodyFace;
104
105 typedef struct ReferenceVert {
106 float pos[3]; /* position relative to com */
107 float mass; /* node mass */
108 } ReferenceVert;
109
110 typedef struct ReferenceState {
111 float com[3]; /* center of mass*/
112 ReferenceVert *ivert; /* list of initial values */
113 } ReferenceState;
114
115 /*private scratch pad for caching and other data only needed when alive*/
116 typedef struct SBScratch {
117 GHash *colliderhash;
118 short needstobuildcollider;
119 short flag;
120 BodyFace *bodyface;
121 int totface;
122 float aabbmin[3], aabbmax[3];
123 ReferenceState Ref;
124 } SBScratch;
125
126 typedef struct SB_thread_context {
127 Scene *scene;
128 Object *ob;
129 float forcetime;
130 float timenow;
131 int ifirst;
132 int ilast;
133 ListBase *effectors;
134 int do_deflector;
135 float fieldfactor;
136 float windfactor;
137 int nr;
138 int tot;
139 } SB_thread_context;
140
141 #define MID_PRESERVE 1
142
143 #define SOFTGOALSNAP 0.999f
144 /* if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp
145 * removes *unnecessary* stiffness from ODE system
146 */
147 #define HEUNWARNLIMIT 1 /* 500 would be fine i think for detecting severe *stiff* stuff */
148
149 #define BSF_INTERSECT 1 /* edge intersects collider face */
150
151 /* private definitions for bodypoint states */
152 #define SBF_DOFUZZY 1 /* Bodypoint do fuzzy */
153 #define SBF_OUTOFCOLLISION 2 /* Bodypoint does not collide */
154
155 #define BFF_INTERSECT 1 /* collider edge intrudes face */
156 #define BFF_CLOSEVERT 2 /* collider vertex repulses face */
157
158 /* humm .. this should be calculated from sb parameters and sizes. */
159 static float SoftHeunTol = 1.0f;
160
161 /* local prototypes */
162 static void free_softbody_intern(SoftBody *sb);
163
164 /*+++ frame based timing +++*/
165
166 /*physical unit of force is [kg * m / sec^2]*/
167
168 /**
169 * Since unit of g is [m/sec^2] and F = mass * g we re-scale unit mass of node to 1 gram
170 * put it to a function here, so we can add user options later without touching simulation code.
171 */
sb_grav_force_scale(Object * UNUSED (ob))172 static float sb_grav_force_scale(Object *UNUSED(ob))
173 {
174 return (0.001f);
175 }
176
177 /**
178 * Re-scaling unit of drag [1 / sec] to somehow reasonable
179 * put it to a function here, so we can add user options later without touching simulation code.
180 */
sb_fric_force_scale(Object * UNUSED (ob))181 static float sb_fric_force_scale(Object *UNUSED(ob))
182 {
183 return (0.01f);
184 }
185
186 /**
187 * Defining the frames to *real* time relation.
188 */
sb_time_scale(Object * ob)189 static float sb_time_scale(Object *ob)
190 {
191 SoftBody *sb = ob->soft; /* is supposed to be there */
192 if (sb) {
193 return (sb->physics_speed);
194 /* hrms .. this could be IPO as well :)
195 * estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)]
196 * 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions] has period 65 frames
197 * theory would give a 50 frames period .. so there must be something inaccurate ..
198 * looking for that (BM). */
199 }
200 return (1.0f);
201 /*
202 * this would be frames/sec independent timing assuming 25 fps is default
203 * but does not work very well with NLA
204 * return (25.0f/scene->r.frs_sec)
205 */
206 }
207 /*--- frame based timing ---*/
208
209 /* helper functions for everything is animatable jow_go_for2_5 +++++++*/
210 /* introducing them here, because i know: steps in properties ( at frame timing )
211 * will cause unwanted responses of the softbody system (which does inter frame calculations )
212 * so first 'cure' would be: interpolate linear in time ..
213 * Q: why do i write this?
214 * A: because it happened once, that some eager coder 'streamlined' code to fail.
215 * We DO linear interpolation for goals .. and i think we should do on animated properties as well
216 */
217
218 /* animate sb->maxgoal, sb->mingoal */
_final_goal(Object * ob,BodyPoint * bp)219 static float _final_goal(Object *ob, BodyPoint *bp) /*jow_go_for2_5 */
220 {
221 float f = -1999.99f;
222 if (ob) {
223 SoftBody *sb = ob->soft; /* is supposed to be there */
224 if (!(ob->softflag & OB_SB_GOAL)) {
225 return (0.0f);
226 }
227 if (sb && bp) {
228 if (bp->goal < 0.0f) {
229 return (0.0f);
230 }
231 f = sb->mingoal + bp->goal * fabsf(sb->maxgoal - sb->mingoal);
232 f = pow(f, 4.0f);
233 return f;
234 }
235 }
236 CLOG_ERROR(&LOG, "sb or bp == NULL");
237 return f; /*using crude but spot able values some times helps debuggin */
238 }
239
_final_mass(Object * ob,BodyPoint * bp)240 static float _final_mass(Object *ob, BodyPoint *bp)
241 {
242 if (ob) {
243 SoftBody *sb = ob->soft; /* is supposed to be there */
244 if (sb && bp) {
245 return (bp->mass * sb->nodemass);
246 }
247 }
248 CLOG_ERROR(&LOG, "sb or bp == NULL");
249 return 1.0f;
250 }
251 /* helper functions for everything is animateble jow_go_for2_5 ------*/
252
253 /*+++ collider caching and dicing +++*/
254
255 /*
256 * for each target object/face the axis aligned bounding box (AABB) is stored
257 * faces parallel to global axes
258 * so only simple "value" in [min, max] checks are used
259 * float operations still
260 */
261
262 /* just an ID here to reduce the prob for killing objects
263 * ob->sumohandle points to we should not kill :)
264 */
265 static const int CCD_SAFETY = 190561;
266
267 typedef struct ccdf_minmax {
268 float minx, miny, minz, maxx, maxy, maxz;
269 } ccdf_minmax;
270
271 typedef struct ccd_Mesh {
272 int mvert_num, tri_num;
273 const MVert *mvert;
274 const MVert *mprevvert;
275 const MVertTri *tri;
276 int safety;
277 ccdf_minmax *mima;
278 /* Axis Aligned Bounding Box AABB */
279 float bbmin[3];
280 float bbmax[3];
281 } ccd_Mesh;
282
ccd_mesh_make(Object * ob)283 static ccd_Mesh *ccd_mesh_make(Object *ob)
284 {
285 CollisionModifierData *cmd;
286 ccd_Mesh *pccd_M = NULL;
287 ccdf_minmax *mima;
288 const MVertTri *vt;
289 float hull;
290 int i;
291
292 cmd = (CollisionModifierData *)BKE_modifiers_findby_type(ob, eModifierType_Collision);
293
294 /* first some paranoia checks */
295 if (!cmd) {
296 return NULL;
297 }
298 if (!cmd->mvert_num || !cmd->tri_num) {
299 return NULL;
300 }
301
302 pccd_M = MEM_mallocN(sizeof(ccd_Mesh), "ccd_Mesh");
303 pccd_M->mvert_num = cmd->mvert_num;
304 pccd_M->tri_num = cmd->tri_num;
305 pccd_M->safety = CCD_SAFETY;
306 pccd_M->bbmin[0] = pccd_M->bbmin[1] = pccd_M->bbmin[2] = 1e30f;
307 pccd_M->bbmax[0] = pccd_M->bbmax[1] = pccd_M->bbmax[2] = -1e30f;
308 pccd_M->mprevvert = NULL;
309
310 /* blow it up with forcefield ranges */
311 hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft);
312
313 /* alloc and copy verts*/
314 pccd_M->mvert = MEM_dupallocN(cmd->xnew);
315 /* note that xnew coords are already in global space, */
316 /* determine the ortho BB */
317 for (i = 0; i < pccd_M->mvert_num; i++) {
318 const float *v;
319
320 /* evaluate limits */
321 v = pccd_M->mvert[i].co;
322 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
323 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
324 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
325
326 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
327 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
328 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
329 }
330 /* alloc and copy faces*/
331 pccd_M->tri = MEM_dupallocN(cmd->tri);
332
333 /* OBBs for idea1 */
334 pccd_M->mima = MEM_mallocN(sizeof(ccdf_minmax) * pccd_M->tri_num, "ccd_Mesh_Faces_mima");
335
336 /* anyhoo we need to walk the list of faces and find OBB they live in */
337 for (i = 0, mima = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
338 const float *v;
339
340 mima->minx = mima->miny = mima->minz = 1e30f;
341 mima->maxx = mima->maxy = mima->maxz = -1e30f;
342
343 v = pccd_M->mvert[vt->tri[0]].co;
344 mima->minx = min_ff(mima->minx, v[0] - hull);
345 mima->miny = min_ff(mima->miny, v[1] - hull);
346 mima->minz = min_ff(mima->minz, v[2] - hull);
347 mima->maxx = max_ff(mima->maxx, v[0] + hull);
348 mima->maxy = max_ff(mima->maxy, v[1] + hull);
349 mima->maxz = max_ff(mima->maxz, v[2] + hull);
350
351 v = pccd_M->mvert[vt->tri[1]].co;
352 mima->minx = min_ff(mima->minx, v[0] - hull);
353 mima->miny = min_ff(mima->miny, v[1] - hull);
354 mima->minz = min_ff(mima->minz, v[2] - hull);
355 mima->maxx = max_ff(mima->maxx, v[0] + hull);
356 mima->maxy = max_ff(mima->maxy, v[1] + hull);
357 mima->maxz = max_ff(mima->maxz, v[2] + hull);
358
359 v = pccd_M->mvert[vt->tri[2]].co;
360 mima->minx = min_ff(mima->minx, v[0] - hull);
361 mima->miny = min_ff(mima->miny, v[1] - hull);
362 mima->minz = min_ff(mima->minz, v[2] - hull);
363 mima->maxx = max_ff(mima->maxx, v[0] + hull);
364 mima->maxy = max_ff(mima->maxy, v[1] + hull);
365 mima->maxz = max_ff(mima->maxz, v[2] + hull);
366 }
367
368 return pccd_M;
369 }
ccd_mesh_update(Object * ob,ccd_Mesh * pccd_M)370 static void ccd_mesh_update(Object *ob, ccd_Mesh *pccd_M)
371 {
372 CollisionModifierData *cmd;
373 ccdf_minmax *mima;
374 const MVertTri *vt;
375 float hull;
376 int i;
377
378 cmd = (CollisionModifierData *)BKE_modifiers_findby_type(ob, eModifierType_Collision);
379
380 /* first some paranoia checks */
381 if (!cmd) {
382 return;
383 }
384 if (!cmd->mvert_num || !cmd->tri_num) {
385 return;
386 }
387
388 if ((pccd_M->mvert_num != cmd->mvert_num) || (pccd_M->tri_num != cmd->tri_num)) {
389 return;
390 }
391
392 pccd_M->bbmin[0] = pccd_M->bbmin[1] = pccd_M->bbmin[2] = 1e30f;
393 pccd_M->bbmax[0] = pccd_M->bbmax[1] = pccd_M->bbmax[2] = -1e30f;
394
395 /* blow it up with forcefield ranges */
396 hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft);
397
398 /* rotate current to previous */
399 if (pccd_M->mprevvert) {
400 MEM_freeN((void *)pccd_M->mprevvert);
401 }
402 pccd_M->mprevvert = pccd_M->mvert;
403 /* alloc and copy verts*/
404 pccd_M->mvert = MEM_dupallocN(cmd->xnew);
405 /* note that xnew coords are already in global space, */
406 /* determine the ortho BB */
407 for (i = 0; i < pccd_M->mvert_num; i++) {
408 const float *v;
409
410 /* evaluate limits */
411 v = pccd_M->mvert[i].co;
412 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
413 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
414 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
415
416 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
417 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
418 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
419
420 /* evaluate limits */
421 v = pccd_M->mprevvert[i].co;
422 pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull);
423 pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull);
424 pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull);
425
426 pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull);
427 pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull);
428 pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull);
429 }
430
431 /* anyhoo we need to walk the list of faces and find OBB they live in */
432 for (i = 0, mima = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
433 const float *v;
434
435 mima->minx = mima->miny = mima->minz = 1e30f;
436 mima->maxx = mima->maxy = mima->maxz = -1e30f;
437
438 /* mvert */
439 v = pccd_M->mvert[vt->tri[0]].co;
440 mima->minx = min_ff(mima->minx, v[0] - hull);
441 mima->miny = min_ff(mima->miny, v[1] - hull);
442 mima->minz = min_ff(mima->minz, v[2] - hull);
443 mima->maxx = max_ff(mima->maxx, v[0] + hull);
444 mima->maxy = max_ff(mima->maxy, v[1] + hull);
445 mima->maxz = max_ff(mima->maxz, v[2] + hull);
446
447 v = pccd_M->mvert[vt->tri[1]].co;
448 mima->minx = min_ff(mima->minx, v[0] - hull);
449 mima->miny = min_ff(mima->miny, v[1] - hull);
450 mima->minz = min_ff(mima->minz, v[2] - hull);
451 mima->maxx = max_ff(mima->maxx, v[0] + hull);
452 mima->maxy = max_ff(mima->maxy, v[1] + hull);
453 mima->maxz = max_ff(mima->maxz, v[2] + hull);
454
455 v = pccd_M->mvert[vt->tri[2]].co;
456 mima->minx = min_ff(mima->minx, v[0] - hull);
457 mima->miny = min_ff(mima->miny, v[1] - hull);
458 mima->minz = min_ff(mima->minz, v[2] - hull);
459 mima->maxx = max_ff(mima->maxx, v[0] + hull);
460 mima->maxy = max_ff(mima->maxy, v[1] + hull);
461 mima->maxz = max_ff(mima->maxz, v[2] + hull);
462
463 /* mprevvert */
464 v = pccd_M->mprevvert[vt->tri[0]].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 v = pccd_M->mprevvert[vt->tri[1]].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 v = pccd_M->mprevvert[vt->tri[2]].co;
481 mima->minx = min_ff(mima->minx, v[0] - hull);
482 mima->miny = min_ff(mima->miny, v[1] - hull);
483 mima->minz = min_ff(mima->minz, v[2] - hull);
484 mima->maxx = max_ff(mima->maxx, v[0] + hull);
485 mima->maxy = max_ff(mima->maxy, v[1] + hull);
486 mima->maxz = max_ff(mima->maxz, v[2] + hull);
487 }
488 }
489
ccd_mesh_free(ccd_Mesh * ccdm)490 static void ccd_mesh_free(ccd_Mesh *ccdm)
491 {
492 /* Make sure we're not nuking objects we don't know. */
493 if (ccdm && (ccdm->safety == CCD_SAFETY)) {
494 MEM_freeN((void *)ccdm->mvert);
495 MEM_freeN((void *)ccdm->tri);
496 if (ccdm->mprevvert) {
497 MEM_freeN((void *)ccdm->mprevvert);
498 }
499 MEM_freeN(ccdm->mima);
500 MEM_freeN(ccdm);
501 ccdm = NULL;
502 }
503 }
504
ccd_build_deflector_hash_single(GHash * hash,Object * ob)505 static void ccd_build_deflector_hash_single(GHash *hash, Object *ob)
506 {
507 /* only with deflecting set */
508 if (ob->pd && ob->pd->deflect) {
509 void **val_p;
510 if (!BLI_ghash_ensure_p(hash, ob, &val_p)) {
511 ccd_Mesh *ccdmesh = ccd_mesh_make(ob);
512 *val_p = ccdmesh;
513 }
514 }
515 }
516
517 /**
518 * \note collection overrides scene when not NULL.
519 */
ccd_build_deflector_hash(Depsgraph * depsgraph,Collection * collection,Object * vertexowner,GHash * hash)520 static void ccd_build_deflector_hash(Depsgraph *depsgraph,
521 Collection *collection,
522 Object *vertexowner,
523 GHash *hash)
524 {
525 if (!hash) {
526 return;
527 }
528
529 unsigned int numobjects;
530 Object **objects = BKE_collision_objects_create(
531 depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
532
533 for (int i = 0; i < numobjects; i++) {
534 Object *ob = objects[i];
535
536 if (ob->type == OB_MESH) {
537 ccd_build_deflector_hash_single(hash, ob);
538 }
539 }
540
541 BKE_collision_objects_free(objects);
542 }
543
ccd_update_deflector_hash_single(GHash * hash,Object * ob)544 static void ccd_update_deflector_hash_single(GHash *hash, Object *ob)
545 {
546 if (ob->pd && ob->pd->deflect) {
547 ccd_Mesh *ccdmesh = BLI_ghash_lookup(hash, ob);
548 if (ccdmesh) {
549 ccd_mesh_update(ob, ccdmesh);
550 }
551 }
552 }
553
554 /**
555 * \note collection overrides scene when not NULL.
556 */
ccd_update_deflector_hash(Depsgraph * depsgraph,Collection * collection,Object * vertexowner,GHash * hash)557 static void ccd_update_deflector_hash(Depsgraph *depsgraph,
558 Collection *collection,
559 Object *vertexowner,
560 GHash *hash)
561 {
562 if ((!hash) || (!vertexowner)) {
563 return;
564 }
565
566 unsigned int numobjects;
567 Object **objects = BKE_collision_objects_create(
568 depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
569
570 for (int i = 0; i < numobjects; i++) {
571 Object *ob = objects[i];
572
573 if (ob->type == OB_MESH) {
574 ccd_update_deflector_hash_single(hash, ob);
575 }
576 }
577
578 BKE_collision_objects_free(objects);
579 }
580
581 /*--- collider caching and dicing ---*/
582
count_mesh_quads(Mesh * me)583 static int count_mesh_quads(Mesh *me)
584 {
585 int a, result = 0;
586 const MPoly *mp = me->mpoly;
587
588 if (mp) {
589 for (a = me->totpoly; a > 0; a--, mp++) {
590 if (mp->totloop == 4) {
591 result++;
592 }
593 }
594 }
595 return result;
596 }
597
add_mesh_quad_diag_springs(Object * ob)598 static void add_mesh_quad_diag_springs(Object *ob)
599 {
600 Mesh *me = ob->data;
601 /*BodyPoint *bp;*/ /*UNUSED*/
602 int a;
603
604 if (ob->soft) {
605 int nofquads;
606 // float s_shear = ob->soft->shearstiff*ob->soft->shearstiff;
607
608 nofquads = count_mesh_quads(me);
609 if (nofquads) {
610 const MLoop *mloop = me->mloop;
611 const MPoly *mp = me->mpoly;
612 BodySpring *bs;
613
614 /* resize spring-array to hold additional quad springs */
615 ob->soft->bspring = MEM_recallocN(ob->soft->bspring,
616 sizeof(BodySpring) * (ob->soft->totspring + nofquads * 2));
617
618 /* fill the tail */
619 a = 0;
620 bs = &ob->soft->bspring[ob->soft->totspring];
621 /*bp= ob->soft->bpoint; */ /*UNUSED*/
622 for (a = me->totpoly; a > 0; a--, mp++) {
623 if (mp->totloop == 4) {
624 bs->v1 = mloop[mp->loopstart + 0].v;
625 bs->v2 = mloop[mp->loopstart + 2].v;
626 bs->springtype = SB_STIFFQUAD;
627 bs++;
628 bs->v1 = mloop[mp->loopstart + 1].v;
629 bs->v2 = mloop[mp->loopstart + 3].v;
630 bs->springtype = SB_STIFFQUAD;
631 bs++;
632 }
633 }
634
635 /* now we can announce new springs */
636 ob->soft->totspring += nofquads * 2;
637 }
638 }
639 }
640
add_2nd_order_roller(Object * ob,float UNUSED (stiffness),int * counter,int addsprings)641 static void add_2nd_order_roller(Object *ob, float UNUSED(stiffness), int *counter, int addsprings)
642 {
643 /*assume we have a softbody*/
644 SoftBody *sb = ob->soft; /* is supposed to be there */
645 BodyPoint *bp, *bpo;
646 BodySpring *bs, *bs2, *bs3 = NULL;
647 int a, b, c, notthis = 0, v0;
648 if (!sb->bspring) {
649 return;
650 } /* we are 2nd order here so 1rst should have been build :) */
651 /* first run counting second run adding */
652 *counter = 0;
653 if (addsprings) {
654 bs3 = ob->soft->bspring + ob->soft->totspring;
655 }
656 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
657 /*scan for neighborhood*/
658 bpo = NULL;
659 v0 = (sb->totpoint - a);
660 for (b = bp->nofsprings; b > 0; b--) {
661 bs = sb->bspring + bp->springs[b - 1];
662 /* Nasty thing here that springs have two ends
663 * so here we have to make sure we examine the other */
664 if (v0 == bs->v1) {
665 bpo = sb->bpoint + bs->v2;
666 notthis = bs->v2;
667 }
668 else {
669 if (v0 == bs->v2) {
670 bpo = sb->bpoint + bs->v1;
671 notthis = bs->v1;
672 }
673 else {
674 CLOG_ERROR(&LOG, "oops we should not get here");
675 }
676 }
677 if (bpo) { /* so now we have a 2nd order humpdidump */
678 for (c = bpo->nofsprings; c > 0; c--) {
679 bs2 = sb->bspring + bpo->springs[c - 1];
680 if ((bs2->v1 != notthis) && (bs2->v1 > v0)) {
681 (*counter)++; /*hit */
682 if (addsprings) {
683 bs3->v1 = v0;
684 bs3->v2 = bs2->v1;
685 bs3->springtype = SB_BEND;
686 bs3++;
687 }
688 }
689 if ((bs2->v2 != notthis) && (bs2->v2 > v0)) {
690 (*counter)++; /* hit */
691 if (addsprings) {
692 bs3->v1 = v0;
693 bs3->v2 = bs2->v2;
694 bs3->springtype = SB_BEND;
695 bs3++;
696 }
697 }
698 }
699 }
700 }
701 /*scan for neighborhood done*/
702 }
703 }
704
add_2nd_order_springs(Object * ob,float stiffness)705 static void add_2nd_order_springs(Object *ob, float stiffness)
706 {
707 int counter = 0;
708 BodySpring *bs_new;
709 stiffness *= stiffness;
710
711 add_2nd_order_roller(ob, stiffness, &counter, 0); /* counting */
712 if (counter) {
713 /* resize spring-array to hold additional springs */
714 bs_new = MEM_callocN((ob->soft->totspring + counter) * sizeof(BodySpring), "bodyspring");
715 memcpy(bs_new, ob->soft->bspring, (ob->soft->totspring) * sizeof(BodySpring));
716
717 if (ob->soft->bspring) {
718 MEM_freeN(ob->soft->bspring);
719 }
720 ob->soft->bspring = bs_new;
721
722 add_2nd_order_roller(ob, stiffness, &counter, 1); /* adding */
723 ob->soft->totspring += counter;
724 }
725 }
726
add_bp_springlist(BodyPoint * bp,int springID)727 static void add_bp_springlist(BodyPoint *bp, int springID)
728 {
729 int *newlist;
730
731 if (bp->springs == NULL) {
732 bp->springs = MEM_callocN(sizeof(int), "bpsprings");
733 bp->springs[0] = springID;
734 bp->nofsprings = 1;
735 }
736 else {
737 bp->nofsprings++;
738 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings");
739 memcpy(newlist, bp->springs, (bp->nofsprings - 1) * sizeof(int));
740 MEM_freeN(bp->springs);
741 bp->springs = newlist;
742 bp->springs[bp->nofsprings - 1] = springID;
743 }
744 }
745
746 /**
747 * Do this once when sb is build it is `O(N^2)`
748 * so scanning for springs every iteration is too expensive.
749 */
build_bps_springlist(Object * ob)750 static void build_bps_springlist(Object *ob)
751 {
752 SoftBody *sb = ob->soft; /* is supposed to be there */
753 BodyPoint *bp;
754 BodySpring *bs;
755 int a, b;
756
757 if (sb == NULL) {
758 return; /* paranoia check */
759 }
760
761 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
762 /* throw away old list */
763 if (bp->springs) {
764 MEM_freeN(bp->springs);
765 bp->springs = NULL;
766 }
767 /* scan for attached inner springs */
768 for (b = sb->totspring, bs = sb->bspring; b > 0; b--, bs++) {
769 if (((sb->totpoint - a) == bs->v1)) {
770 add_bp_springlist(bp, sb->totspring - b);
771 }
772 if (((sb->totpoint - a) == bs->v2)) {
773 add_bp_springlist(bp, sb->totspring - b);
774 }
775 } /*for springs*/
776 } /*for bp*/
777 }
778
calculate_collision_balls(Object * ob)779 static void calculate_collision_balls(Object *ob)
780 {
781 SoftBody *sb = ob->soft; /* is supposed to be there */
782 BodyPoint *bp;
783 BodySpring *bs;
784 int a, b, akku_count;
785 float min, max, akku;
786
787 if (sb == NULL) {
788 return; /* paranoia check */
789 }
790
791 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
792 bp->colball = 0;
793 akku = 0.0f;
794 akku_count = 0;
795 min = 1e22f;
796 max = -1e22f;
797 /* first estimation based on attached */
798 for (b = bp->nofsprings; b > 0; b--) {
799 bs = sb->bspring + bp->springs[b - 1];
800 if (bs->springtype == SB_EDGE) {
801 akku += bs->len;
802 akku_count++;
803 min = min_ff(bs->len, min);
804 max = max_ff(bs->len, max);
805 }
806 }
807
808 if (akku_count > 0) {
809 if (sb->sbc_mode == SBC_MODE_MANUAL) {
810 bp->colball = sb->colball;
811 }
812 if (sb->sbc_mode == SBC_MODE_AVG) {
813 bp->colball = akku / (float)akku_count * sb->colball;
814 }
815 if (sb->sbc_mode == SBC_MODE_MIN) {
816 bp->colball = min * sb->colball;
817 }
818 if (sb->sbc_mode == SBC_MODE_MAX) {
819 bp->colball = max * sb->colball;
820 }
821 if (sb->sbc_mode == SBC_MODE_AVGMINMAX) {
822 bp->colball = (min + max) / 2.0f * sb->colball;
823 }
824 }
825 else {
826 bp->colball = 0;
827 }
828 } /*for bp*/
829 }
830
831 /* creates new softbody if didn't exist yet, makes new points and springs arrays */
renew_softbody(Scene * scene,Object * ob,int totpoint,int totspring)832 static void renew_softbody(Scene *scene, Object *ob, int totpoint, int totspring)
833 {
834 SoftBody *sb;
835 int i;
836 short softflag;
837 if (ob->soft == NULL) {
838 ob->soft = sbNew(scene);
839 }
840 else {
841 free_softbody_intern(ob->soft);
842 }
843 sb = ob->soft;
844 softflag = ob->softflag;
845
846 if (totpoint) {
847 sb->totpoint = totpoint;
848 sb->totspring = totspring;
849
850 sb->bpoint = MEM_mallocN(totpoint * sizeof(BodyPoint), "bodypoint");
851 if (totspring) {
852 sb->bspring = MEM_mallocN(totspring * sizeof(BodySpring), "bodyspring");
853 }
854
855 /* initialize BodyPoint array */
856 for (i = 0; i < totpoint; i++) {
857 BodyPoint *bp = &sb->bpoint[i];
858
859 /* hum as far as i see this is overridden by _final_goal() now jow_go_for2_5 */
860 /* sadly breaks compatibility with older versions */
861 /* but makes goals behave the same for meshes, lattices and curves */
862 if (softflag & OB_SB_GOAL) {
863 bp->goal = sb->defgoal;
864 }
865 else {
866 bp->goal = 0.0f;
867 /* so this will definily be below SOFTGOALSNAP */
868 }
869
870 bp->nofsprings = 0;
871 bp->springs = NULL;
872 bp->choke = 0.0f;
873 bp->choke2 = 0.0f;
874 bp->frozen = 1.0f;
875 bp->colball = 0.0f;
876 bp->loc_flag = 0;
877 bp->springweight = 1.0f;
878 bp->mass = 1.0f;
879 }
880 }
881 }
882
free_softbody_baked(SoftBody * sb)883 static void free_softbody_baked(SoftBody *sb)
884 {
885 SBVertex *key;
886 int k;
887
888 for (k = 0; k < sb->totkey; k++) {
889 key = *(sb->keys + k);
890 if (key) {
891 MEM_freeN(key);
892 }
893 }
894 if (sb->keys) {
895 MEM_freeN(sb->keys);
896 }
897
898 sb->keys = NULL;
899 sb->totkey = 0;
900 }
free_scratch(SoftBody * sb)901 static void free_scratch(SoftBody *sb)
902 {
903 if (sb->scratch) {
904 /* todo make sure everything is cleaned up nicly */
905 if (sb->scratch->colliderhash) {
906 BLI_ghash_free(sb->scratch->colliderhash,
907 NULL,
908 (GHashValFreeFP)ccd_mesh_free); /*this hoepfully will free all caches*/
909 sb->scratch->colliderhash = NULL;
910 }
911 if (sb->scratch->bodyface) {
912 MEM_freeN(sb->scratch->bodyface);
913 }
914 if (sb->scratch->Ref.ivert) {
915 MEM_freeN(sb->scratch->Ref.ivert);
916 }
917 MEM_freeN(sb->scratch);
918 sb->scratch = NULL;
919 }
920 }
921
922 /* only frees internal data */
free_softbody_intern(SoftBody * sb)923 static void free_softbody_intern(SoftBody *sb)
924 {
925 if (sb) {
926 int a;
927 BodyPoint *bp;
928
929 if (sb->bpoint) {
930 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
931 /* free spring list */
932 if (bp->springs != NULL) {
933 MEM_freeN(bp->springs);
934 }
935 }
936 MEM_freeN(sb->bpoint);
937 }
938
939 if (sb->bspring) {
940 MEM_freeN(sb->bspring);
941 }
942
943 sb->totpoint = sb->totspring = 0;
944 sb->bpoint = NULL;
945 sb->bspring = NULL;
946
947 free_scratch(sb);
948 free_softbody_baked(sb);
949 }
950 }
951
952 /* ************ dynamics ********** */
953
954 /* the most general (micro physics correct) way to do collision
955 * (only needs the current particle position)
956 *
957 * it actually checks if the particle intrudes a short range force field generated
958 * by the faces of the target object and returns a force to drive the particle out
959 * the strength of the field grows exponentially if the particle is on the 'wrong' side of the face
960 * 'wrong' side : projection to the face normal is negative (all referred to a vertex in the face)
961 *
962 * flaw of this: 'fast' particles as well as 'fast' colliding faces
963 * give a 'tunnel' effect such that the particle passes through the force field
964 * without ever 'seeing' it
965 * this is fully compliant to heisenberg: h >= fuzzy(location) * fuzzy(time)
966 * besides our h is way larger than in QM because forces propagate way slower here
967 * we have to deal with fuzzy(time) in the range of 1/25 seconds (typical frame rate)
968 * yup collision targets are not known here any better
969 * and 1/25 second is very long compared to real collision events
970 * Q: why not use 'simple' collision here like bouncing back a particle
971 * --> reverting is velocity on the face normal
972 * A: because our particles are not alone here
973 * and need to tell their neighbors exactly what happens via spring forces
974 * unless sbObjectStep( .. ) is called on sub frame timing level
975 * BTW that also questions the use of a 'implicit' solvers on softbodies
976 * since that would only valid for 'slow' moving collision targets and dito particles
977 */
978
979 /* +++ dependency information functions*/
980
981 /**
982 * \note collection overrides scene when not NULL.
983 */
query_external_colliders(Depsgraph * depsgraph,Collection * collection)984 static int query_external_colliders(Depsgraph *depsgraph, Collection *collection)
985 {
986 unsigned int numobjects;
987 Object **objects = BKE_collision_objects_create(
988 depsgraph, NULL, collection, &numobjects, eModifierType_Collision);
989 BKE_collision_objects_free(objects);
990
991 return (numobjects != 0);
992 }
993 /* --- dependency information functions*/
994
995 /* +++ the aabb "force" section*/
sb_detect_aabb_collisionCached(float UNUSED (force[3]),struct Object * vertexowner,float UNUSED (time))996 static int sb_detect_aabb_collisionCached(float UNUSED(force[3]),
997 struct Object *vertexowner,
998 float UNUSED(time))
999 {
1000 Object *ob;
1001 SoftBody *sb = vertexowner->soft;
1002 GHash *hash;
1003 GHashIterator *ihash;
1004 float aabbmin[3], aabbmax[3];
1005 int deflected = 0;
1006 #if 0
1007 int a;
1008 #endif
1009
1010 if ((sb == NULL) || (sb->scratch == NULL)) {
1011 return 0;
1012 }
1013 copy_v3_v3(aabbmin, sb->scratch->aabbmin);
1014 copy_v3_v3(aabbmax, sb->scratch->aabbmax);
1015
1016 hash = vertexowner->soft->scratch->colliderhash;
1017 ihash = BLI_ghashIterator_new(hash);
1018 while (!BLI_ghashIterator_done(ihash)) {
1019
1020 ccd_Mesh *ccdm = BLI_ghashIterator_getValue(ihash);
1021 ob = BLI_ghashIterator_getKey(ihash);
1022 {
1023 /* only with deflecting set */
1024 if (ob->pd && ob->pd->deflect) {
1025 if (ccdm) {
1026 if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) ||
1027 (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) ||
1028 (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2])) {
1029 /* boxes don't intersect */
1030 BLI_ghashIterator_step(ihash);
1031 continue;
1032 }
1033
1034 /* so now we have the 2 boxes overlapping */
1035 /* forces actually not used */
1036 deflected = 2;
1037 }
1038 else {
1039 /*aye that should be cached*/
1040 CLOG_ERROR(&LOG, "missing cache error");
1041 BLI_ghashIterator_step(ihash);
1042 continue;
1043 }
1044 } /* if (ob->pd && ob->pd->deflect) */
1045 BLI_ghashIterator_step(ihash);
1046 }
1047 } /* while () */
1048 BLI_ghashIterator_free(ihash);
1049 return deflected;
1050 }
1051 /* --- the aabb section*/
1052
1053 /* +++ the face external section*/
sb_detect_face_pointCached(const float face_v1[3],const float face_v2[3],const float face_v3[3],float * damp,float force[3],struct Object * vertexowner,float time)1054 static int sb_detect_face_pointCached(const float face_v1[3],
1055 const float face_v2[3],
1056 const float face_v3[3],
1057 float *damp,
1058 float force[3],
1059 struct Object *vertexowner,
1060 float time)
1061 {
1062 Object *ob;
1063 GHash *hash;
1064 GHashIterator *ihash;
1065 float nv1[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1066 float facedist, outerfacethickness, tune = 10.f;
1067 int a, deflected = 0;
1068
1069 aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]);
1070 aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]);
1071 aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]);
1072 aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]);
1073 aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]);
1074 aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]);
1075
1076 /* calculate face normal once again SIGH */
1077 sub_v3_v3v3(edge1, face_v1, face_v2);
1078 sub_v3_v3v3(edge2, face_v3, face_v2);
1079 cross_v3_v3v3(d_nvect, edge2, edge1);
1080 normalize_v3(d_nvect);
1081
1082 hash = vertexowner->soft->scratch->colliderhash;
1083 ihash = BLI_ghashIterator_new(hash);
1084 while (!BLI_ghashIterator_done(ihash)) {
1085
1086 ccd_Mesh *ccdm = BLI_ghashIterator_getValue(ihash);
1087 ob = BLI_ghashIterator_getKey(ihash);
1088 {
1089 /* only with deflecting set */
1090 if (ob->pd && ob->pd->deflect) {
1091 const MVert *mvert = NULL;
1092 const MVert *mprevvert = NULL;
1093 if (ccdm) {
1094 mvert = ccdm->mvert;
1095 a = ccdm->mvert_num;
1096 mprevvert = ccdm->mprevvert;
1097 outerfacethickness = ob->pd->pdef_sboft;
1098 if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) ||
1099 (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) ||
1100 (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2])) {
1101 /* boxes don't intersect */
1102 BLI_ghashIterator_step(ihash);
1103 continue;
1104 }
1105 }
1106 else {
1107 /*aye that should be cached*/
1108 CLOG_ERROR(&LOG, "missing cache error");
1109 BLI_ghashIterator_step(ihash);
1110 continue;
1111 }
1112
1113 /* use mesh*/
1114 if (mvert) {
1115 while (a) {
1116 copy_v3_v3(nv1, mvert[a - 1].co);
1117 if (mprevvert) {
1118 mul_v3_fl(nv1, time);
1119 madd_v3_v3fl(nv1, mprevvert[a - 1].co, 1.0f - time);
1120 }
1121 /* origin to face_v2*/
1122 sub_v3_v3(nv1, face_v2);
1123 facedist = dot_v3v3(nv1, d_nvect);
1124 if (fabsf(facedist) < outerfacethickness) {
1125 if (isect_point_tri_prism_v3(nv1, face_v1, face_v2, face_v3)) {
1126 float df;
1127 if (facedist > 0) {
1128 df = (outerfacethickness - facedist) / outerfacethickness;
1129 }
1130 else {
1131 df = (outerfacethickness + facedist) / outerfacethickness;
1132 }
1133
1134 *damp = df * tune * ob->pd->pdef_sbdamp;
1135
1136 df = 0.01f * expf(-100.0f * df);
1137 madd_v3_v3fl(force, d_nvect, -df);
1138 deflected = 3;
1139 }
1140 }
1141 a--;
1142 } /* while (a)*/
1143 } /* if (mvert) */
1144 } /* if (ob->pd && ob->pd->deflect) */
1145 BLI_ghashIterator_step(ihash);
1146 }
1147 } /* while () */
1148 BLI_ghashIterator_free(ihash);
1149 return deflected;
1150 }
1151
sb_detect_face_collisionCached(const float face_v1[3],const float face_v2[3],const float face_v3[3],float * damp,float force[3],struct Object * vertexowner,float time)1152 static int sb_detect_face_collisionCached(const float face_v1[3],
1153 const float face_v2[3],
1154 const float face_v3[3],
1155 float *damp,
1156 float force[3],
1157 struct Object *vertexowner,
1158 float time)
1159 {
1160 Object *ob;
1161 GHash *hash;
1162 GHashIterator *ihash;
1163 float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1164 float t, tune = 10.0f;
1165 int a, deflected = 0;
1166
1167 aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]);
1168 aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]);
1169 aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]);
1170 aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]);
1171 aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]);
1172 aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]);
1173
1174 hash = vertexowner->soft->scratch->colliderhash;
1175 ihash = BLI_ghashIterator_new(hash);
1176 while (!BLI_ghashIterator_done(ihash)) {
1177
1178 ccd_Mesh *ccdm = BLI_ghashIterator_getValue(ihash);
1179 ob = BLI_ghashIterator_getKey(ihash);
1180 {
1181 /* only with deflecting set */
1182 if (ob->pd && ob->pd->deflect) {
1183 const MVert *mvert = NULL;
1184 const MVert *mprevvert = NULL;
1185 const MVertTri *vt = NULL;
1186 const ccdf_minmax *mima = NULL;
1187
1188 if (ccdm) {
1189 mvert = ccdm->mvert;
1190 vt = ccdm->tri;
1191 mprevvert = ccdm->mprevvert;
1192 mima = ccdm->mima;
1193 a = ccdm->tri_num;
1194
1195 if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) ||
1196 (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) ||
1197 (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2])) {
1198 /* boxes don't intersect */
1199 BLI_ghashIterator_step(ihash);
1200 continue;
1201 }
1202 }
1203 else {
1204 /*aye that should be cached*/
1205 CLOG_ERROR(&LOG, "missing cache error");
1206 BLI_ghashIterator_step(ihash);
1207 continue;
1208 }
1209
1210 /* use mesh*/
1211 while (a--) {
1212 if ((aabbmax[0] < mima->minx) || (aabbmin[0] > mima->maxx) ||
1213 (aabbmax[1] < mima->miny) || (aabbmin[1] > mima->maxy) ||
1214 (aabbmax[2] < mima->minz) || (aabbmin[2] > mima->maxz)) {
1215 mima++;
1216 vt++;
1217 continue;
1218 }
1219
1220 if (mvert) {
1221
1222 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1223 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1224 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1225
1226 if (mprevvert) {
1227 mul_v3_fl(nv1, time);
1228 madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1229
1230 mul_v3_fl(nv2, time);
1231 madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1232
1233 mul_v3_fl(nv3, time);
1234 madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1235 }
1236 }
1237
1238 /* switch origin to be nv2*/
1239 sub_v3_v3v3(edge1, nv1, nv2);
1240 sub_v3_v3v3(edge2, nv3, nv2);
1241 cross_v3_v3v3(d_nvect, edge2, edge1);
1242 normalize_v3(d_nvect);
1243 if (isect_line_segment_tri_v3(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) ||
1244 isect_line_segment_tri_v3(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) ||
1245 isect_line_segment_tri_v3(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL)) {
1246 madd_v3_v3fl(force, d_nvect, -0.5f);
1247 *damp = tune * ob->pd->pdef_sbdamp;
1248 deflected = 2;
1249 }
1250 mima++;
1251 vt++;
1252 } /* while a */
1253 } /* if (ob->pd && ob->pd->deflect) */
1254 BLI_ghashIterator_step(ihash);
1255 }
1256 } /* while () */
1257 BLI_ghashIterator_free(ihash);
1258 return deflected;
1259 }
1260
scan_for_ext_face_forces(Object * ob,float timenow)1261 static void scan_for_ext_face_forces(Object *ob, float timenow)
1262 {
1263 SoftBody *sb = ob->soft;
1264 BodyFace *bf;
1265 int a;
1266 float damp = 0.0f, choke = 1.0f;
1267 float tune = -10.0f;
1268 float feedback[3];
1269
1270 if (sb && sb->scratch->totface) {
1271
1272 bf = sb->scratch->bodyface;
1273 for (a = 0; a < sb->scratch->totface; a++, bf++) {
1274 bf->ext_force[0] = bf->ext_force[1] = bf->ext_force[2] = 0.0f;
1275 /*+++edges intruding*/
1276 bf->flag &= ~BFF_INTERSECT;
1277 zero_v3(feedback);
1278 if (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos,
1279 sb->bpoint[bf->v2].pos,
1280 sb->bpoint[bf->v3].pos,
1281 &damp,
1282 feedback,
1283 ob,
1284 timenow)) {
1285 madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
1286 madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
1287 madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
1288 // madd_v3_v3fl(bf->ext_force, feedback, tune);
1289 bf->flag |= BFF_INTERSECT;
1290 choke = min_ff(max_ff(damp, choke), 1.0f);
1291 }
1292 /*---edges intruding*/
1293
1294 /*+++ close vertices*/
1295 if ((bf->flag & BFF_INTERSECT) == 0) {
1296 bf->flag &= ~BFF_CLOSEVERT;
1297 tune = -1.0f;
1298 zero_v3(feedback);
1299 if (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos,
1300 sb->bpoint[bf->v2].pos,
1301 sb->bpoint[bf->v3].pos,
1302 &damp,
1303 feedback,
1304 ob,
1305 timenow)) {
1306 madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
1307 madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
1308 madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
1309 // madd_v3_v3fl(bf->ext_force, feedback, tune);
1310 bf->flag |= BFF_CLOSEVERT;
1311 choke = min_ff(max_ff(damp, choke), 1.0f);
1312 }
1313 }
1314 /*--- close vertices*/
1315 }
1316 bf = sb->scratch->bodyface;
1317 for (a = 0; a < sb->scratch->totface; a++, bf++) {
1318 if ((bf->flag & BFF_INTERSECT) || (bf->flag & BFF_CLOSEVERT)) {
1319 sb->bpoint[bf->v1].choke2 = max_ff(sb->bpoint[bf->v1].choke2, choke);
1320 sb->bpoint[bf->v2].choke2 = max_ff(sb->bpoint[bf->v2].choke2, choke);
1321 sb->bpoint[bf->v3].choke2 = max_ff(sb->bpoint[bf->v3].choke2, choke);
1322 }
1323 }
1324 }
1325 }
1326
1327 /* --- the face external section*/
1328
1329 /* +++ the spring external section*/
1330
sb_detect_edge_collisionCached(const float edge_v1[3],const float edge_v2[3],float * damp,float force[3],struct Object * vertexowner,float time)1331 static int sb_detect_edge_collisionCached(const float edge_v1[3],
1332 const float edge_v2[3],
1333 float *damp,
1334 float force[3],
1335 struct Object *vertexowner,
1336 float time)
1337 {
1338 Object *ob;
1339 GHash *hash;
1340 GHashIterator *ihash;
1341 float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3];
1342 float t, el;
1343 int a, deflected = 0;
1344
1345 minmax_v3v3_v3(aabbmin, aabbmax, edge_v1);
1346 minmax_v3v3_v3(aabbmin, aabbmax, edge_v2);
1347
1348 el = len_v3v3(edge_v1, edge_v2);
1349
1350 hash = vertexowner->soft->scratch->colliderhash;
1351 ihash = BLI_ghashIterator_new(hash);
1352 while (!BLI_ghashIterator_done(ihash)) {
1353
1354 ccd_Mesh *ccdm = BLI_ghashIterator_getValue(ihash);
1355 ob = BLI_ghashIterator_getKey(ihash);
1356 {
1357 /* only with deflecting set */
1358 if (ob->pd && ob->pd->deflect) {
1359 const MVert *mvert = NULL;
1360 const MVert *mprevvert = NULL;
1361 const MVertTri *vt = NULL;
1362 const ccdf_minmax *mima = NULL;
1363
1364 if (ccdm) {
1365 mvert = ccdm->mvert;
1366 mprevvert = ccdm->mprevvert;
1367 vt = ccdm->tri;
1368 mima = ccdm->mima;
1369 a = ccdm->tri_num;
1370
1371 if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) ||
1372 (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) ||
1373 (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2])) {
1374 /* boxes don't intersect */
1375 BLI_ghashIterator_step(ihash);
1376 continue;
1377 }
1378 }
1379 else {
1380 /*aye that should be cached*/
1381 CLOG_ERROR(&LOG, "missing cache error");
1382 BLI_ghashIterator_step(ihash);
1383 continue;
1384 }
1385
1386 /* use mesh*/
1387 while (a--) {
1388 if ((aabbmax[0] < mima->minx) || (aabbmin[0] > mima->maxx) ||
1389 (aabbmax[1] < mima->miny) || (aabbmin[1] > mima->maxy) ||
1390 (aabbmax[2] < mima->minz) || (aabbmin[2] > mima->maxz)) {
1391 mima++;
1392 vt++;
1393 continue;
1394 }
1395
1396 if (mvert) {
1397
1398 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1399 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1400 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1401
1402 if (mprevvert) {
1403 mul_v3_fl(nv1, time);
1404 madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1405
1406 mul_v3_fl(nv2, time);
1407 madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1408
1409 mul_v3_fl(nv3, time);
1410 madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1411 }
1412 }
1413
1414 /* switch origin to be nv2*/
1415 sub_v3_v3v3(edge1, nv1, nv2);
1416 sub_v3_v3v3(edge2, nv3, nv2);
1417
1418 cross_v3_v3v3(d_nvect, edge2, edge1);
1419 normalize_v3(d_nvect);
1420 if (isect_line_segment_tri_v3(edge_v1, edge_v2, nv1, nv2, nv3, &t, NULL)) {
1421 float v1[3], v2[3];
1422 float intrusiondepth, i1, i2;
1423 sub_v3_v3v3(v1, edge_v1, nv2);
1424 sub_v3_v3v3(v2, edge_v2, nv2);
1425 i1 = dot_v3v3(v1, d_nvect);
1426 i2 = dot_v3v3(v2, d_nvect);
1427 intrusiondepth = -min_ff(i1, i2) / el;
1428 madd_v3_v3fl(force, d_nvect, intrusiondepth);
1429 *damp = ob->pd->pdef_sbdamp;
1430 deflected = 2;
1431 }
1432
1433 mima++;
1434 vt++;
1435 } /* while a */
1436 } /* if (ob->pd && ob->pd->deflect) */
1437 BLI_ghashIterator_step(ihash);
1438 }
1439 } /* while () */
1440 BLI_ghashIterator_free(ihash);
1441 return deflected;
1442 }
1443
_scan_for_ext_spring_forces(Scene * scene,Object * ob,float timenow,int ifirst,int ilast,struct ListBase * effectors)1444 static void _scan_for_ext_spring_forces(
1445 Scene *scene, Object *ob, float timenow, int ifirst, int ilast, struct ListBase *effectors)
1446 {
1447 SoftBody *sb = ob->soft;
1448 int a;
1449 float damp;
1450 float feedback[3];
1451
1452 if (sb && sb->totspring) {
1453 for (a = ifirst; a < ilast; a++) {
1454 BodySpring *bs = &sb->bspring[a];
1455 bs->ext_force[0] = bs->ext_force[1] = bs->ext_force[2] = 0.0f;
1456 feedback[0] = feedback[1] = feedback[2] = 0.0f;
1457 bs->flag &= ~BSF_INTERSECT;
1458
1459 if (bs->springtype == SB_EDGE) {
1460 /* +++ springs colliding */
1461 if (ob->softflag & OB_SB_EDGECOLL) {
1462 if (sb_detect_edge_collisionCached(
1463 sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos, &damp, feedback, ob, timenow)) {
1464 add_v3_v3(bs->ext_force, feedback);
1465 bs->flag |= BSF_INTERSECT;
1466 // bs->cf=damp;
1467 bs->cf = sb->choke * 0.01f;
1468 }
1469 }
1470 /* ---- springs colliding */
1471
1472 /* +++ springs seeing wind ... n stuff depending on their orientation*/
1473 /* note we don't use sb->mediafrict but use sb->aeroedge for magnitude of effect*/
1474 if (sb->aeroedge) {
1475 float vel[3], sp[3], pr[3], force[3];
1476 float f, windfactor = 0.25f;
1477 /*see if we have wind*/
1478 if (effectors) {
1479 EffectedPoint epoint;
1480 float speed[3] = {0.0f, 0.0f, 0.0f};
1481 float pos[3];
1482 mid_v3_v3v3(pos, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos);
1483 mid_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec);
1484 pd_point_from_soft(scene, pos, vel, -1, &epoint);
1485 BKE_effectors_apply(
1486 effectors, NULL, sb->effector_weights, &epoint, force, NULL, speed);
1487
1488 mul_v3_fl(speed, windfactor);
1489 add_v3_v3(vel, speed);
1490 }
1491 /* media in rest */
1492 else {
1493 add_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec);
1494 }
1495 f = normalize_v3(vel);
1496 f = -0.0001f * f * f * sb->aeroedge;
1497 /* (todo) add a nice angle dependent function done for now BUT */
1498 /* still there could be some nice drag/lift function, but who needs it */
1499
1500 sub_v3_v3v3(sp, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos);
1501 project_v3_v3v3(pr, vel, sp);
1502 sub_v3_v3(vel, pr);
1503 normalize_v3(vel);
1504 if (ob->softflag & OB_SB_AERO_ANGLE) {
1505 normalize_v3(sp);
1506 madd_v3_v3fl(bs->ext_force, vel, f * (1.0f - fabsf(dot_v3v3(vel, sp))));
1507 }
1508 else {
1509 madd_v3_v3fl(bs->ext_force, vel, f); /* to keep compatible with 2.45 release files */
1510 }
1511 }
1512 /* --- springs seeing wind */
1513 }
1514 }
1515 }
1516 }
1517
exec_scan_for_ext_spring_forces(void * data)1518 static void *exec_scan_for_ext_spring_forces(void *data)
1519 {
1520 SB_thread_context *pctx = (SB_thread_context *)data;
1521 _scan_for_ext_spring_forces(
1522 pctx->scene, pctx->ob, pctx->timenow, pctx->ifirst, pctx->ilast, pctx->effectors);
1523 return NULL;
1524 }
1525
sb_sfesf_threads_run(struct Depsgraph * depsgraph,Scene * scene,struct Object * ob,float timenow,int totsprings,int * UNUSED (ptr_to_break_func (void)))1526 static void sb_sfesf_threads_run(struct Depsgraph *depsgraph,
1527 Scene *scene,
1528 struct Object *ob,
1529 float timenow,
1530 int totsprings,
1531 int *UNUSED(ptr_to_break_func(void)))
1532 {
1533 ListBase threads;
1534 SB_thread_context *sb_threads;
1535 int i, totthread, left, dec;
1536
1537 /* wild guess .. may increase with better thread management 'above'
1538 * or even be UI option sb->spawn_cf_threads_nopts */
1539 int lowsprings = 100;
1540
1541 ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, ob->soft->effector_weights);
1542
1543 /* figure the number of threads while preventing pretty pointless threading overhead */
1544 totthread = BKE_scene_num_threads(scene);
1545 /* what if we got zillions of CPUs running but less to spread*/
1546 while ((totsprings / totthread < lowsprings) && (totthread > 1)) {
1547 totthread--;
1548 }
1549
1550 sb_threads = MEM_callocN(sizeof(SB_thread_context) * totthread, "SBSpringsThread");
1551 memset(sb_threads, 0, sizeof(SB_thread_context) * totthread);
1552 left = totsprings;
1553 dec = totsprings / totthread + 1;
1554 for (i = 0; i < totthread; i++) {
1555 sb_threads[i].scene = scene;
1556 sb_threads[i].ob = ob;
1557 sb_threads[i].forcetime = 0.0; /* not used here */
1558 sb_threads[i].timenow = timenow;
1559 sb_threads[i].ilast = left;
1560 left = left - dec;
1561 if (left > 0) {
1562 sb_threads[i].ifirst = left;
1563 }
1564 else {
1565 sb_threads[i].ifirst = 0;
1566 }
1567 sb_threads[i].effectors = effectors;
1568 sb_threads[i].do_deflector = false; /* not used here */
1569 sb_threads[i].fieldfactor = 0.0f; /* not used here */
1570 sb_threads[i].windfactor = 0.0f; /* not used here */
1571 sb_threads[i].nr = i;
1572 sb_threads[i].tot = totthread;
1573 }
1574 if (totthread > 1) {
1575 BLI_threadpool_init(&threads, exec_scan_for_ext_spring_forces, totthread);
1576
1577 for (i = 0; i < totthread; i++) {
1578 BLI_threadpool_insert(&threads, &sb_threads[i]);
1579 }
1580
1581 BLI_threadpool_end(&threads);
1582 }
1583 else {
1584 exec_scan_for_ext_spring_forces(&sb_threads[0]);
1585 }
1586 /* clean up */
1587 MEM_freeN(sb_threads);
1588
1589 BKE_effectors_free(effectors);
1590 }
1591
1592 /* --- the spring external section*/
1593
choose_winner(float * w,float * pos,float * a,float * b,float * c,float * ca,float * cb,float * cc)1594 static int choose_winner(
1595 float *w, float *pos, float *a, float *b, float *c, float *ca, float *cb, float *cc)
1596 {
1597 float mindist, cp;
1598 int winner = 1;
1599 mindist = fabsf(dot_v3v3(pos, a));
1600
1601 cp = fabsf(dot_v3v3(pos, b));
1602 if (mindist < cp) {
1603 mindist = cp;
1604 winner = 2;
1605 }
1606
1607 cp = fabsf(dot_v3v3(pos, c));
1608 if (mindist < cp) {
1609 mindist = cp;
1610 winner = 3;
1611 }
1612 switch (winner) {
1613 case 1:
1614 copy_v3_v3(w, ca);
1615 break;
1616 case 2:
1617 copy_v3_v3(w, cb);
1618 break;
1619 case 3:
1620 copy_v3_v3(w, cc);
1621 }
1622 return winner;
1623 }
1624
sb_detect_vertex_collisionCached(float opco[3],float facenormal[3],float * damp,float force[3],struct Object * vertexowner,float time,float vel[3],float * intrusion)1625 static int sb_detect_vertex_collisionCached(float opco[3],
1626 float facenormal[3],
1627 float *damp,
1628 float force[3],
1629 struct Object *vertexowner,
1630 float time,
1631 float vel[3],
1632 float *intrusion)
1633 {
1634 Object *ob = NULL;
1635 GHash *hash;
1636 GHashIterator *ihash;
1637 float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], dv1[3], ve[3],
1638 avel[3] = {0.0, 0.0, 0.0}, vv1[3], vv2[3], vv3[3], coledge[3] = {0.0f, 0.0f, 0.0f},
1639 mindistedge = 1000.0f, outerforceaccu[3], innerforceaccu[3], facedist,
1640 /* n_mag, */ /* UNUSED */ force_mag_norm, minx, miny, minz, maxx, maxy, maxz,
1641 innerfacethickness = -0.5f, outerfacethickness = 0.2f, ee = 5.0f, ff = 0.1f, fa = 1;
1642 int a, deflected = 0, cavel = 0, ci = 0;
1643 /* init */
1644 *intrusion = 0.0f;
1645 hash = vertexowner->soft->scratch->colliderhash;
1646 ihash = BLI_ghashIterator_new(hash);
1647 outerforceaccu[0] = outerforceaccu[1] = outerforceaccu[2] = 0.0f;
1648 innerforceaccu[0] = innerforceaccu[1] = innerforceaccu[2] = 0.0f;
1649 /* go */
1650 while (!BLI_ghashIterator_done(ihash)) {
1651
1652 ccd_Mesh *ccdm = BLI_ghashIterator_getValue(ihash);
1653 ob = BLI_ghashIterator_getKey(ihash);
1654 {
1655 /* only with deflecting set */
1656 if (ob->pd && ob->pd->deflect) {
1657 const MVert *mvert = NULL;
1658 const MVert *mprevvert = NULL;
1659 const MVertTri *vt = NULL;
1660 const ccdf_minmax *mima = NULL;
1661
1662 if (ccdm) {
1663 mvert = ccdm->mvert;
1664 mprevvert = ccdm->mprevvert;
1665 vt = ccdm->tri;
1666 mima = ccdm->mima;
1667 a = ccdm->tri_num;
1668
1669 minx = ccdm->bbmin[0];
1670 miny = ccdm->bbmin[1];
1671 minz = ccdm->bbmin[2];
1672
1673 maxx = ccdm->bbmax[0];
1674 maxy = ccdm->bbmax[1];
1675 maxz = ccdm->bbmax[2];
1676
1677 if ((opco[0] < minx) || (opco[1] < miny) || (opco[2] < minz) || (opco[0] > maxx) ||
1678 (opco[1] > maxy) || (opco[2] > maxz)) {
1679 /* outside the padded boundbox --> collision object is too far away */
1680 BLI_ghashIterator_step(ihash);
1681 continue;
1682 }
1683 }
1684 else {
1685 /*aye that should be cached*/
1686 CLOG_ERROR(&LOG, "missing cache error");
1687 BLI_ghashIterator_step(ihash);
1688 continue;
1689 }
1690
1691 /* do object level stuff */
1692 /* need to have user control for that since it depends on model scale */
1693 innerfacethickness = -ob->pd->pdef_sbift;
1694 outerfacethickness = ob->pd->pdef_sboft;
1695 fa = (ff * outerfacethickness - outerfacethickness);
1696 fa *= fa;
1697 fa = 1.0f / fa;
1698 avel[0] = avel[1] = avel[2] = 0.0f;
1699 /* use mesh*/
1700 while (a--) {
1701 if ((opco[0] < mima->minx) || (opco[0] > mima->maxx) || (opco[1] < mima->miny) ||
1702 (opco[1] > mima->maxy) || (opco[2] < mima->minz) || (opco[2] > mima->maxz)) {
1703 mima++;
1704 vt++;
1705 continue;
1706 }
1707
1708 if (mvert) {
1709
1710 copy_v3_v3(nv1, mvert[vt->tri[0]].co);
1711 copy_v3_v3(nv2, mvert[vt->tri[1]].co);
1712 copy_v3_v3(nv3, mvert[vt->tri[2]].co);
1713
1714 if (mprevvert) {
1715 /* Grab the average speed of the collider vertices before we spoil nvX
1716 * humm could be done once a SB steps but then we' need to store that too
1717 * since the AABB reduced probability to get here drastically
1718 * it might be a nice tradeoff CPU <--> memory.
1719 */
1720 sub_v3_v3v3(vv1, nv1, mprevvert[vt->tri[0]].co);
1721 sub_v3_v3v3(vv2, nv2, mprevvert[vt->tri[1]].co);
1722 sub_v3_v3v3(vv3, nv3, mprevvert[vt->tri[2]].co);
1723
1724 mul_v3_fl(nv1, time);
1725 madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
1726
1727 mul_v3_fl(nv2, time);
1728 madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
1729
1730 mul_v3_fl(nv3, time);
1731 madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
1732 }
1733 }
1734
1735 /* switch origin to be nv2*/
1736 sub_v3_v3v3(edge1, nv1, nv2);
1737 sub_v3_v3v3(edge2, nv3, nv2);
1738 /* Abuse dv1 to have vertex in question at *origin* of triangle. */
1739 sub_v3_v3v3(dv1, opco, nv2);
1740
1741 cross_v3_v3v3(d_nvect, edge2, edge1);
1742 /* n_mag = */ /* UNUSED */ normalize_v3(d_nvect);
1743 facedist = dot_v3v3(dv1, d_nvect);
1744 /* so rules are */
1745
1746 if ((facedist > innerfacethickness) && (facedist < outerfacethickness)) {
1747 if (isect_point_tri_prism_v3(opco, nv1, nv2, nv3)) {
1748 force_mag_norm = (float)exp(-ee * facedist);
1749 if (facedist > outerfacethickness * ff) {
1750 force_mag_norm = (float)force_mag_norm * fa * (facedist - outerfacethickness) *
1751 (facedist - outerfacethickness);
1752 }
1753 *damp = ob->pd->pdef_sbdamp;
1754 if (facedist > 0.0f) {
1755 *damp *= (1.0f - facedist / outerfacethickness);
1756 madd_v3_v3fl(outerforceaccu, d_nvect, force_mag_norm);
1757 deflected = 3;
1758 }
1759 else {
1760 madd_v3_v3fl(innerforceaccu, d_nvect, force_mag_norm);
1761 if (deflected < 2) {
1762 deflected = 2;
1763 }
1764 }
1765 if ((mprevvert) && (*damp > 0.0f)) {
1766 choose_winner(ve, opco, nv1, nv2, nv3, vv1, vv2, vv3);
1767 add_v3_v3(avel, ve);
1768 cavel++;
1769 }
1770 *intrusion += facedist;
1771 ci++;
1772 }
1773 }
1774
1775 mima++;
1776 vt++;
1777 } /* while a */
1778 } /* if (ob->pd && ob->pd->deflect) */
1779 BLI_ghashIterator_step(ihash);
1780 }
1781 } /* while () */
1782
1783 if (deflected == 1) { /* no face but 'outer' edge cylinder sees vert */
1784 force_mag_norm = (float)exp(-ee * mindistedge);
1785 if (mindistedge > outerfacethickness * ff) {
1786 force_mag_norm = (float)force_mag_norm * fa * (mindistedge - outerfacethickness) *
1787 (mindistedge - outerfacethickness);
1788 }
1789 madd_v3_v3fl(force, coledge, force_mag_norm);
1790 *damp = ob->pd->pdef_sbdamp;
1791 if (mindistedge > 0.0f) {
1792 *damp *= (1.0f - mindistedge / outerfacethickness);
1793 }
1794 }
1795 if (deflected == 2) { /* face inner detected */
1796 add_v3_v3(force, innerforceaccu);
1797 }
1798 if (deflected == 3) { /* face outer detected */
1799 add_v3_v3(force, outerforceaccu);
1800 }
1801
1802 BLI_ghashIterator_free(ihash);
1803 if (cavel) {
1804 mul_v3_fl(avel, 1.0f / (float)cavel);
1805 }
1806 copy_v3_v3(vel, avel);
1807 if (ci) {
1808 *intrusion /= ci;
1809 }
1810 if (deflected) {
1811 normalize_v3_v3(facenormal, force);
1812 }
1813 return deflected;
1814 }
1815
1816 /* sandbox to plug in various deflection algos */
sb_deflect_face(Object * ob,float * actpos,float * facenormal,float * force,float * cf,float time,float * vel,float * intrusion)1817 static int sb_deflect_face(Object *ob,
1818 float *actpos,
1819 float *facenormal,
1820 float *force,
1821 float *cf,
1822 float time,
1823 float *vel,
1824 float *intrusion)
1825 {
1826 float s_actpos[3];
1827 int deflected;
1828 copy_v3_v3(s_actpos, actpos);
1829 deflected = sb_detect_vertex_collisionCached(
1830 s_actpos, facenormal, cf, force, ob, time, vel, intrusion);
1831 #if 0
1832 deflected = sb_detect_vertex_collisionCachedEx(
1833 s_actpos, facenormal, cf, force, ob, time, vel, intrusion);
1834 #endif
1835 return deflected;
1836 }
1837
1838 /* hiding this for now .. but the jacobian may pop up on other tasks .. so i'd like to keep it */
1839 #if 0
1840 static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len, float factor)
1841 {
1842 float m, delta_ij;
1843 int i, j;
1844 if (L < len) {
1845 for (i = 0; i < 3; i++) {
1846 for (j = 0; j < 3; j++) {
1847 delta_ij = (i == j ? (1.0f) : (0.0f));
1848 m = factor * (dir[i] * dir[j] + (1 - L / len) * (delta_ij - dir[i] * dir[j]));
1849 EIG_linear_solver_matrix_add(ia + i, op + ic + j, m);
1850 }
1851 }
1852 }
1853 else {
1854 for (i = 0; i < 3; i++) {
1855 for (j = 0; j < 3; j++) {
1856 m = factor * dir[i] * dir[j];
1857 EIG_linear_solver_matrix_add(ia + i, op + ic + j, m);
1858 }
1859 }
1860 }
1861 }
1862
1863 static void dfdx_goal(int ia, int ic, int op, float factor)
1864 {
1865 int i;
1866 for (i = 0; i < 3; i++) {
1867 EIG_linear_solver_matrix_add(ia + i, op + ic + i, factor);
1868 }
1869 }
1870
1871 static void dfdv_goal(int ia, int ic, float factor)
1872 {
1873 int i;
1874 for (i = 0; i < 3; i++) {
1875 EIG_linear_solver_matrix_add(ia + i, ic + i, factor);
1876 }
1877 }
1878 #endif /* if 0 */
1879
sb_spring_force(Object * ob,int bpi,BodySpring * bs,float iks,float UNUSED (forcetime))1880 static void sb_spring_force(
1881 Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime))
1882 {
1883 SoftBody *sb = ob->soft; /* is supposed to be there */
1884 BodyPoint *bp1, *bp2;
1885
1886 float dir[3], dvel[3];
1887 float distance, forcefactor, kd, absvel, projvel, kw;
1888 #if 0 /* UNUSED */
1889 int ia, ic;
1890 #endif
1891 /* prepare depending on which side of the spring we are on */
1892 if (bpi == bs->v1) {
1893 bp1 = &sb->bpoint[bs->v1];
1894 bp2 = &sb->bpoint[bs->v2];
1895 #if 0 /* UNUSED */
1896 ia = 3 * bs->v1;
1897 ic = 3 * bs->v2;
1898 #endif
1899 }
1900 else if (bpi == bs->v2) {
1901 bp1 = &sb->bpoint[bs->v2];
1902 bp2 = &sb->bpoint[bs->v1];
1903 #if 0 /* UNUSED */
1904 ia = 3 * bs->v2;
1905 ic = 3 * bs->v1;
1906 #endif
1907 }
1908 else {
1909 /* TODO make this debug option */
1910 CLOG_WARN(&LOG, "bodypoint <bpi> is not attached to spring <*bs>");
1911 return;
1912 }
1913
1914 /* do bp1 <--> bp2 elastic */
1915 sub_v3_v3v3(dir, bp1->pos, bp2->pos);
1916 distance = normalize_v3(dir);
1917 if (bs->len < distance) {
1918 iks = 1.0f / (1.0f - sb->inspring) - 1.0f; /* inner spring constants function */
1919 }
1920 else {
1921 iks = 1.0f / (1.0f - sb->inpush) - 1.0f; /* inner spring constants function */
1922 }
1923
1924 if (bs->len > 0.0f) { /* check for degenerated springs */
1925 forcefactor = iks / bs->len;
1926 }
1927 else {
1928 forcefactor = iks;
1929 }
1930 kw = (bp1->springweight + bp2->springweight) / 2.0f;
1931 kw = kw * kw;
1932 kw = kw * kw;
1933 switch (bs->springtype) {
1934 case SB_EDGE:
1935 case SB_HANDLE:
1936 forcefactor *= kw;
1937 break;
1938 case SB_BEND:
1939 forcefactor *= sb->secondspring * kw;
1940 break;
1941 case SB_STIFFQUAD:
1942 forcefactor *= sb->shearstiff * sb->shearstiff * kw;
1943 break;
1944 default:
1945 break;
1946 }
1947
1948 madd_v3_v3fl(bp1->force, dir, (bs->len - distance) * forcefactor);
1949
1950 /* do bp1 <--> bp2 viscous */
1951 sub_v3_v3v3(dvel, bp1->vec, bp2->vec);
1952 kd = sb->infrict * sb_fric_force_scale(ob);
1953 absvel = normalize_v3(dvel);
1954 projvel = dot_v3v3(dir, dvel);
1955 kd *= absvel * projvel;
1956 madd_v3_v3fl(bp1->force, dir, -kd);
1957 }
1958
1959 /* since this is definitely the most CPU consuming task here .. try to spread it */
1960 /* core function _softbody_calc_forces_slice_in_a_thread */
1961 /* result is int to be able to flag user break */
_softbody_calc_forces_slice_in_a_thread(Scene * scene,Object * ob,float forcetime,float timenow,int ifirst,int ilast,int * UNUSED (ptr_to_break_func (void)),ListBase * effectors,int do_deflector,float fieldfactor,float windfactor)1962 static int _softbody_calc_forces_slice_in_a_thread(Scene *scene,
1963 Object *ob,
1964 float forcetime,
1965 float timenow,
1966 int ifirst,
1967 int ilast,
1968 int *UNUSED(ptr_to_break_func(void)),
1969 ListBase *effectors,
1970 int do_deflector,
1971 float fieldfactor,
1972 float windfactor)
1973 {
1974 float iks;
1975 int bb, do_selfcollision, do_springcollision, do_aero;
1976 int number_of_points_here = ilast - ifirst;
1977 SoftBody *sb = ob->soft; /* is supposed to be there */
1978 BodyPoint *bp;
1979
1980 /* initialize */
1981 if (sb) {
1982 /* check conditions for various options */
1983 /* +++ could be done on object level to squeeze out the last bits of it */
1984 do_selfcollision = ((ob->softflag & OB_SB_EDGES) && (sb->bspring) &&
1985 (ob->softflag & OB_SB_SELF));
1986 do_springcollision = do_deflector && (ob->softflag & OB_SB_EDGES) &&
1987 (ob->softflag & OB_SB_EDGECOLL);
1988 do_aero = ((sb->aeroedge) && (ob->softflag & OB_SB_EDGES));
1989 /* --- could be done on object level to squeeze out the last bits of it */
1990 }
1991 else {
1992 CLOG_ERROR(&LOG, "expected a SB here");
1993 return 999;
1994 }
1995
1996 /* debugerin */
1997 if (sb->totpoint < ifirst) {
1998 printf("Aye 998");
1999 return 998;
2000 }
2001 /* debugerin */
2002
2003 bp = &sb->bpoint[ifirst];
2004 for (bb = number_of_points_here; bb > 0; bb--, bp++) {
2005 /* clear forces accumulator */
2006 bp->force[0] = bp->force[1] = bp->force[2] = 0.0;
2007 /* naive ball self collision */
2008 /* needs to be done if goal snaps or not */
2009 if (do_selfcollision) {
2010 int attached;
2011 BodyPoint *obp;
2012 BodySpring *bs;
2013 int c, b;
2014 float velcenter[3], dvel[3], def[3];
2015 float distance;
2016 float compare;
2017 float bstune = sb->ballstiff;
2018
2019 /* Running in a slice we must not assume anything done with obp
2020 * neither alter the data of obp. */
2021 for (c = sb->totpoint, obp = sb->bpoint; c > 0; c--, obp++) {
2022 compare = (obp->colball + bp->colball);
2023 sub_v3_v3v3(def, bp->pos, obp->pos);
2024 /* rather check the AABBoxes before ever calculating the real distance */
2025 /* mathematically it is completely nuts, but performance is pretty much (3) times faster */
2026 if ((fabsf(def[0]) > compare) || (fabsf(def[1]) > compare) || (fabsf(def[2]) > compare)) {
2027 continue;
2028 }
2029 distance = normalize_v3(def);
2030 if (distance < compare) {
2031 /* exclude body points attached with a spring */
2032 attached = 0;
2033 for (b = obp->nofsprings; b > 0; b--) {
2034 bs = sb->bspring + obp->springs[b - 1];
2035 if ((ilast - bb == bs->v2) || (ilast - bb == bs->v1)) {
2036 attached = 1;
2037 continue;
2038 }
2039 }
2040 if (!attached) {
2041 float f = bstune / (distance) + bstune / (compare * compare) * distance -
2042 2.0f * bstune / compare;
2043
2044 mid_v3_v3v3(velcenter, bp->vec, obp->vec);
2045 sub_v3_v3v3(dvel, velcenter, bp->vec);
2046 mul_v3_fl(dvel, _final_mass(ob, bp));
2047
2048 madd_v3_v3fl(bp->force, def, f * (1.0f - sb->balldamp));
2049 madd_v3_v3fl(bp->force, dvel, sb->balldamp);
2050 }
2051 }
2052 }
2053 }
2054 /* naive ball self collision done */
2055
2056 if (_final_goal(ob, bp) < SOFTGOALSNAP) { /* omit this bp when it snaps */
2057 float auxvect[3];
2058 float velgoal[3];
2059
2060 /* do goal stuff */
2061 if (ob->softflag & OB_SB_GOAL) {
2062 /* true elastic goal */
2063 float ks, kd;
2064 sub_v3_v3v3(auxvect, bp->pos, bp->origT);
2065 ks = 1.0f / (1.0f - _final_goal(ob, bp) * sb->goalspring) - 1.0f;
2066 bp->force[0] += -ks * (auxvect[0]);
2067 bp->force[1] += -ks * (auxvect[1]);
2068 bp->force[2] += -ks * (auxvect[2]);
2069
2070 /* calculate damping forces generated by goals*/
2071 sub_v3_v3v3(velgoal, bp->origS, bp->origE);
2072 kd = sb->goalfrict * sb_fric_force_scale(ob);
2073 add_v3_v3v3(auxvect, velgoal, bp->vec);
2074
2075 if (forcetime >
2076 0.0f) { /* make sure friction does not become rocket motor on time reversal */
2077 bp->force[0] -= kd * (auxvect[0]);
2078 bp->force[1] -= kd * (auxvect[1]);
2079 bp->force[2] -= kd * (auxvect[2]);
2080 }
2081 else {
2082 bp->force[0] -= kd * (velgoal[0] - bp->vec[0]);
2083 bp->force[1] -= kd * (velgoal[1] - bp->vec[1]);
2084 bp->force[2] -= kd * (velgoal[2] - bp->vec[2]);
2085 }
2086 }
2087 /* done goal stuff */
2088
2089 /* gravitation */
2090 if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
2091 float gravity[3];
2092 copy_v3_v3(gravity, scene->physics_settings.gravity);
2093
2094 /* Individual mass of node here. */
2095 mul_v3_fl(gravity,
2096 sb_grav_force_scale(ob) * _final_mass(ob, bp) *
2097 sb->effector_weights->global_gravity);
2098
2099 add_v3_v3(bp->force, gravity);
2100 }
2101
2102 /* particle field & vortex */
2103 if (effectors) {
2104 EffectedPoint epoint;
2105 float kd;
2106 float force[3] = {0.0f, 0.0f, 0.0f};
2107 float speed[3] = {0.0f, 0.0f, 0.0f};
2108
2109 /* just for calling function once */
2110 float eval_sb_fric_force_scale = sb_fric_force_scale(ob);
2111
2112 pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint - bp, &epoint);
2113 BKE_effectors_apply(effectors, NULL, sb->effector_weights, &epoint, force, NULL, speed);
2114
2115 /* apply forcefield*/
2116 mul_v3_fl(force, fieldfactor * eval_sb_fric_force_scale);
2117 add_v3_v3(bp->force, force);
2118
2119 /* BP friction in moving media */
2120 kd = sb->mediafrict * eval_sb_fric_force_scale;
2121 bp->force[0] -= kd * (bp->vec[0] + windfactor * speed[0] / eval_sb_fric_force_scale);
2122 bp->force[1] -= kd * (bp->vec[1] + windfactor * speed[1] / eval_sb_fric_force_scale);
2123 bp->force[2] -= kd * (bp->vec[2] + windfactor * speed[2] / eval_sb_fric_force_scale);
2124 /* now we'll have nice centrifugal effect for vortex */
2125 }
2126 else {
2127 /* BP friction in media (not) moving*/
2128 float kd = sb->mediafrict * sb_fric_force_scale(ob);
2129 /* assume it to be proportional to actual velocity */
2130 bp->force[0] -= bp->vec[0] * kd;
2131 bp->force[1] -= bp->vec[1] * kd;
2132 bp->force[2] -= bp->vec[2] * kd;
2133 /* friction in media done */
2134 }
2135 /* +++cached collision targets */
2136 bp->choke = 0.0f;
2137 bp->choke2 = 0.0f;
2138 bp->loc_flag &= ~SBF_DOFUZZY;
2139 if (do_deflector && !(bp->loc_flag & SBF_OUTOFCOLLISION)) {
2140 float cfforce[3], defforce[3] = {0.0f, 0.0f, 0.0f}, vel[3] = {0.0f, 0.0f, 0.0f},
2141 facenormal[3], cf = 1.0f, intrusion;
2142 float kd = 1.0f;
2143
2144 if (sb_deflect_face(ob, bp->pos, facenormal, defforce, &cf, timenow, vel, &intrusion)) {
2145 if (intrusion < 0.0f) {
2146 sb->scratch->flag |= SBF_DOFUZZY;
2147 bp->loc_flag |= SBF_DOFUZZY;
2148 bp->choke = sb->choke * 0.01f;
2149 }
2150
2151 sub_v3_v3v3(cfforce, bp->vec, vel);
2152 madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
2153
2154 madd_v3_v3fl(bp->force, defforce, kd);
2155 }
2156 }
2157 /* ---cached collision targets */
2158
2159 /* +++springs */
2160 iks = 1.0f / (1.0f - sb->inspring) - 1.0f; /* inner spring constants function */
2161 if (ob->softflag & OB_SB_EDGES) {
2162 if (sb->bspring) { /* spring list exists at all ? */
2163 int b;
2164 BodySpring *bs;
2165 for (b = bp->nofsprings; b > 0; b--) {
2166 bs = sb->bspring + bp->springs[b - 1];
2167 if (do_springcollision || do_aero) {
2168 add_v3_v3(bp->force, bs->ext_force);
2169 if (bs->flag & BSF_INTERSECT) {
2170 bp->choke = bs->cf;
2171 }
2172 }
2173 // sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float forcetime)
2174 sb_spring_force(ob, ilast - bb, bs, iks, forcetime);
2175 } /* loop springs */
2176 } /* existing spring list */
2177 } /*any edges*/
2178 /* ---springs */
2179 } /*omit on snap */
2180 } /*loop all bp's*/
2181 return 0; /*done fine*/
2182 }
2183
exec_softbody_calc_forces(void * data)2184 static void *exec_softbody_calc_forces(void *data)
2185 {
2186 SB_thread_context *pctx = (SB_thread_context *)data;
2187 _softbody_calc_forces_slice_in_a_thread(pctx->scene,
2188 pctx->ob,
2189 pctx->forcetime,
2190 pctx->timenow,
2191 pctx->ifirst,
2192 pctx->ilast,
2193 NULL,
2194 pctx->effectors,
2195 pctx->do_deflector,
2196 pctx->fieldfactor,
2197 pctx->windfactor);
2198 return NULL;
2199 }
2200
sb_cf_threads_run(Scene * scene,Object * ob,float forcetime,float timenow,int totpoint,int * UNUSED (ptr_to_break_func (void)),struct ListBase * effectors,int do_deflector,float fieldfactor,float windfactor)2201 static void sb_cf_threads_run(Scene *scene,
2202 Object *ob,
2203 float forcetime,
2204 float timenow,
2205 int totpoint,
2206 int *UNUSED(ptr_to_break_func(void)),
2207 struct ListBase *effectors,
2208 int do_deflector,
2209 float fieldfactor,
2210 float windfactor)
2211 {
2212 ListBase threads;
2213 SB_thread_context *sb_threads;
2214 int i, totthread, left, dec;
2215
2216 /* wild guess .. may increase with better thread management 'above'
2217 * or even be UI option sb->spawn_cf_threads_nopts. */
2218 int lowpoints = 100;
2219
2220 /* figure the number of threads while preventing pretty pointless threading overhead */
2221 totthread = BKE_scene_num_threads(scene);
2222 /* what if we got zillions of CPUs running but less to spread*/
2223 while ((totpoint / totthread < lowpoints) && (totthread > 1)) {
2224 totthread--;
2225 }
2226
2227 /* printf("sb_cf_threads_run spawning %d threads\n", totthread); */
2228
2229 sb_threads = MEM_callocN(sizeof(SB_thread_context) * totthread, "SBThread");
2230 memset(sb_threads, 0, sizeof(SB_thread_context) * totthread);
2231 left = totpoint;
2232 dec = totpoint / totthread + 1;
2233 for (i = 0; i < totthread; i++) {
2234 sb_threads[i].scene = scene;
2235 sb_threads[i].ob = ob;
2236 sb_threads[i].forcetime = forcetime;
2237 sb_threads[i].timenow = timenow;
2238 sb_threads[i].ilast = left;
2239 left = left - dec;
2240 if (left > 0) {
2241 sb_threads[i].ifirst = left;
2242 }
2243 else {
2244 sb_threads[i].ifirst = 0;
2245 }
2246 sb_threads[i].effectors = effectors;
2247 sb_threads[i].do_deflector = do_deflector;
2248 sb_threads[i].fieldfactor = fieldfactor;
2249 sb_threads[i].windfactor = windfactor;
2250 sb_threads[i].nr = i;
2251 sb_threads[i].tot = totthread;
2252 }
2253
2254 if (totthread > 1) {
2255 BLI_threadpool_init(&threads, exec_softbody_calc_forces, totthread);
2256
2257 for (i = 0; i < totthread; i++) {
2258 BLI_threadpool_insert(&threads, &sb_threads[i]);
2259 }
2260
2261 BLI_threadpool_end(&threads);
2262 }
2263 else {
2264 exec_softbody_calc_forces(&sb_threads[0]);
2265 }
2266 /* clean up */
2267 MEM_freeN(sb_threads);
2268 }
2269
softbody_calc_forces(struct Depsgraph * depsgraph,Scene * scene,Object * ob,float forcetime,float timenow)2270 static void softbody_calc_forces(
2271 struct Depsgraph *depsgraph, Scene *scene, Object *ob, float forcetime, float timenow)
2272 {
2273 /* rule we never alter free variables :bp->vec bp->pos in here !
2274 * this will ruin adaptive stepsize AKA heun! (BM)
2275 */
2276 SoftBody *sb = ob->soft; /* is supposed to be there */
2277 /*BodyPoint *bproot;*/ /* UNUSED */
2278 /* float gravity; */ /* UNUSED */
2279 /* float iks; */
2280 float fieldfactor = -1.0f, windfactor = 0.25;
2281 int do_deflector /*, do_selfcollision*/, do_springcollision, do_aero;
2282
2283 /* gravity = sb->grav * sb_grav_force_scale(ob); */ /* UNUSED */
2284
2285 /* check conditions for various options */
2286 do_deflector = query_external_colliders(depsgraph, sb->collision_group);
2287 #if 0
2288 do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF));
2289 #endif
2290 do_springcollision = do_deflector && (ob->softflag & OB_SB_EDGES) &&
2291 (ob->softflag & OB_SB_EDGECOLL);
2292 do_aero = ((sb->aeroedge) && (ob->softflag & OB_SB_EDGES));
2293
2294 /* iks = 1.0f/(1.0f-sb->inspring)-1.0f; */ /* inner spring constants function */ /* UNUSED */
2295 /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */
2296
2297 if (do_springcollision || do_aero) {
2298 sb_sfesf_threads_run(depsgraph, scene, ob, timenow, sb->totspring, NULL);
2299 }
2300
2301 /* after spring scan because it uses Effoctors too */
2302 ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, sb->effector_weights);
2303
2304 if (do_deflector) {
2305 float defforce[3];
2306 do_deflector = sb_detect_aabb_collisionCached(defforce, ob, timenow);
2307 }
2308
2309 sb_cf_threads_run(scene,
2310 ob,
2311 forcetime,
2312 timenow,
2313 sb->totpoint,
2314 NULL,
2315 effectors,
2316 do_deflector,
2317 fieldfactor,
2318 windfactor);
2319
2320 /* finally add forces caused by face collision */
2321 if (ob->softflag & OB_SB_FACECOLL) {
2322 scan_for_ext_face_forces(ob, timenow);
2323 }
2324
2325 /* finish matrix and solve */
2326 BKE_effectors_free(effectors);
2327 }
2328
softbody_apply_forces(Object * ob,float forcetime,int mode,float * err,int mid_flags)2329 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags)
2330 {
2331 /* time evolution */
2332 /* actually does an explicit euler step mode == 0 */
2333 /* or heun ~ 2nd order runge-kutta steps, mode 1, 2 */
2334 SoftBody *sb = ob->soft; /* is supposed to be there */
2335 BodyPoint *bp;
2336 float dx[3] = {0}, dv[3], aabbmin[3], aabbmax[3], cm[3] = {0.0f, 0.0f, 0.0f};
2337 float timeovermass /*, freezeloc=0.00001f, freezeforce=0.00000000001f*/;
2338 float maxerrpos = 0.0f, maxerrvel = 0.0f;
2339 int a, fuzzy = 0;
2340
2341 forcetime *= sb_time_scale(ob);
2342
2343 aabbmin[0] = aabbmin[1] = aabbmin[2] = 1e20f;
2344 aabbmax[0] = aabbmax[1] = aabbmax[2] = -1e20f;
2345
2346 /* old one with homogeneous masses */
2347 /* claim a minimum mass for vertex */
2348 #if 0
2349 if (sb->nodemass > 0.009999f) {
2350 timeovermass = forcetime / sb->nodemass;
2351 }
2352 else {
2353 timeovermass = forcetime / 0.009999f;
2354 }
2355 #endif
2356
2357 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2358 /* Now we have individual masses. */
2359 /* claim a minimum mass for vertex */
2360 if (_final_mass(ob, bp) > 0.009999f) {
2361 timeovermass = forcetime / _final_mass(ob, bp);
2362 }
2363 else {
2364 timeovermass = forcetime / 0.009999f;
2365 }
2366
2367 if (_final_goal(ob, bp) < SOFTGOALSNAP) {
2368 /* this makes t~ = t */
2369 if (mid_flags & MID_PRESERVE) {
2370 copy_v3_v3(dx, bp->vec);
2371 }
2372
2373 /**
2374 * So here is:
2375 * <pre>
2376 * (v)' = a(cceleration) =
2377 * sum(F_springs)/m + gravitation + some friction forces + more forces.
2378 * </pre>
2379 *
2380 * The ( ... )' operator denotes derivate respective time.
2381 *
2382 * The euler step for velocity then becomes:
2383 * <pre>
2384 * v(t + dt) = v(t) + a(t) * dt
2385 * </pre>
2386 */
2387 mul_v3_fl(bp->force, timeovermass); /* individual mass of node here */
2388 /* some nasty if's to have heun in here too */
2389 copy_v3_v3(dv, bp->force);
2390
2391 if (mode == 1) {
2392 copy_v3_v3(bp->prevvec, bp->vec);
2393 copy_v3_v3(bp->prevdv, dv);
2394 }
2395
2396 if (mode == 2) {
2397 /* be optimistic and execute step */
2398 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]);
2399 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]);
2400 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]);
2401 /* compare euler to heun to estimate error for step sizing */
2402 maxerrvel = max_ff(maxerrvel, fabsf(dv[0] - bp->prevdv[0]));
2403 maxerrvel = max_ff(maxerrvel, fabsf(dv[1] - bp->prevdv[1]));
2404 maxerrvel = max_ff(maxerrvel, fabsf(dv[2] - bp->prevdv[2]));
2405 }
2406 else {
2407 add_v3_v3(bp->vec, bp->force);
2408 }
2409
2410 /* this makes t~ = t+dt */
2411 if (!(mid_flags & MID_PRESERVE)) {
2412 copy_v3_v3(dx, bp->vec);
2413 }
2414
2415 /* so here is (x)'= v(elocity) */
2416 /* the euler step for location then becomes */
2417 /* x(t + dt) = x(t) + v(t~) * dt */
2418 mul_v3_fl(dx, forcetime);
2419
2420 /* the freezer coming sooner or later */
2421 #if 0
2422 if ((dot_v3v3(dx, dx) < freezeloc) && (dot_v3v3(bp->force, bp->force) < freezeforce)) {
2423 bp->frozen /= 2;
2424 }
2425 else {
2426 bp->frozen = min_ff(bp->frozen * 1.05f, 1.0f);
2427 }
2428 mul_v3_fl(dx, bp->frozen);
2429 #endif
2430 /* again some nasty if's to have heun in here too */
2431 if (mode == 1) {
2432 copy_v3_v3(bp->prevpos, bp->pos);
2433 copy_v3_v3(bp->prevdx, dx);
2434 }
2435
2436 if (mode == 2) {
2437 bp->pos[0] = bp->prevpos[0] + 0.5f * (dx[0] + bp->prevdx[0]);
2438 bp->pos[1] = bp->prevpos[1] + 0.5f * (dx[1] + bp->prevdx[1]);
2439 bp->pos[2] = bp->prevpos[2] + 0.5f * (dx[2] + bp->prevdx[2]);
2440 maxerrpos = max_ff(maxerrpos, fabsf(dx[0] - bp->prevdx[0]));
2441 maxerrpos = max_ff(maxerrpos, fabsf(dx[1] - bp->prevdx[1]));
2442 maxerrpos = max_ff(maxerrpos, fabsf(dx[2] - bp->prevdx[2]));
2443
2444 /* bp->choke is set when we need to pull a vertex or edge out of the collider.
2445 * the collider object signals to get out by pushing hard. on the other hand
2446 * we don't want to end up in deep space so we add some <viscosity>
2447 * to balance that out */
2448 if (bp->choke2 > 0.0f) {
2449 mul_v3_fl(bp->vec, (1.0f - bp->choke2));
2450 }
2451 if (bp->choke > 0.0f) {
2452 mul_v3_fl(bp->vec, (1.0f - bp->choke));
2453 }
2454 }
2455 else {
2456 add_v3_v3(bp->pos, dx);
2457 }
2458 } /*snap*/
2459 /* so while we are looping BPs anyway do statistics on the fly */
2460 minmax_v3v3_v3(aabbmin, aabbmax, bp->pos);
2461 if (bp->loc_flag & SBF_DOFUZZY) {
2462 fuzzy = 1;
2463 }
2464 } /*for*/
2465
2466 if (sb->totpoint) {
2467 mul_v3_fl(cm, 1.0f / sb->totpoint);
2468 }
2469 if (sb->scratch) {
2470 copy_v3_v3(sb->scratch->aabbmin, aabbmin);
2471 copy_v3_v3(sb->scratch->aabbmax, aabbmax);
2472 }
2473
2474 if (err) { /* so step size will be controlled by biggest difference in slope */
2475 if (sb->solverflags & SBSO_OLDERR) {
2476 *err = max_ff(maxerrpos, maxerrvel);
2477 }
2478 else {
2479 *err = maxerrpos;
2480 }
2481 // printf("EP %f EV %f\n", maxerrpos, maxerrvel);
2482 if (fuzzy) {
2483 *err /= sb->fuzzyness;
2484 }
2485 }
2486 }
2487
2488 /* used by heun when it overshoots */
softbody_restore_prev_step(Object * ob)2489 static void softbody_restore_prev_step(Object *ob)
2490 {
2491 SoftBody *sb = ob->soft; /* is supposed to be there*/
2492 BodyPoint *bp;
2493 int a;
2494
2495 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2496 copy_v3_v3(bp->vec, bp->prevvec);
2497 copy_v3_v3(bp->pos, bp->prevpos);
2498 }
2499 }
2500
2501 #if 0
2502 static void softbody_store_step(Object *ob)
2503 {
2504 SoftBody *sb = ob->soft; /* is supposed to be there*/
2505 BodyPoint *bp;
2506 int a;
2507
2508 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2509 copy_v3_v3(bp->prevvec, bp->vec);
2510 copy_v3_v3(bp->prevpos, bp->pos);
2511 }
2512 }
2513
2514 /* used by predictors and correctors */
2515 static void softbody_store_state(Object *ob, float *ppos, float *pvel)
2516 {
2517 SoftBody *sb = ob->soft; /* is supposed to be there*/
2518 BodyPoint *bp;
2519 int a;
2520 float *pp = ppos, *pv = pvel;
2521
2522 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2523
2524 copy_v3_v3(pv, bp->vec);
2525 pv += 3;
2526
2527 copy_v3_v3(pp, bp->pos);
2528 pp += 3;
2529 }
2530 }
2531
2532 /* used by predictors and correctors */
2533 static void softbody_retrieve_state(Object *ob, float *ppos, float *pvel)
2534 {
2535 SoftBody *sb = ob->soft; /* is supposed to be there*/
2536 BodyPoint *bp;
2537 int a;
2538 float *pp = ppos, *pv = pvel;
2539
2540 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2541
2542 copy_v3_v3(bp->vec, pv);
2543 pv += 3;
2544
2545 copy_v3_v3(bp->pos, pp);
2546 pp += 3;
2547 }
2548 }
2549
2550 /* used by predictors and correctors */
2551 static void softbody_swap_state(Object *ob, float *ppos, float *pvel)
2552 {
2553 SoftBody *sb = ob->soft; /* is supposed to be there*/
2554 BodyPoint *bp;
2555 int a;
2556 float *pp = ppos, *pv = pvel;
2557 float temp[3];
2558
2559 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2560
2561 copy_v3_v3(temp, bp->vec);
2562 copy_v3_v3(bp->vec, pv);
2563 copy_v3_v3(pv, temp);
2564 pv += 3;
2565
2566 copy_v3_v3(temp, bp->pos);
2567 copy_v3_v3(bp->pos, pp);
2568 copy_v3_v3(pp, temp);
2569 pp += 3;
2570 }
2571 }
2572 #endif
2573
2574 /* care for bodypoints taken out of the 'ordinary' solver step
2575 * because they are screwed to goal by bolts
2576 * they just need to move along with the goal in time
2577 * we need to adjust them on sub frame timing in solver
2578 * so now when frame is done .. put 'em to the position at the end of frame
2579 */
softbody_apply_goalsnap(Object * ob)2580 static void softbody_apply_goalsnap(Object *ob)
2581 {
2582 SoftBody *sb = ob->soft; /* is supposed to be there */
2583 BodyPoint *bp;
2584 int a;
2585
2586 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2587 if (_final_goal(ob, bp) >= SOFTGOALSNAP) {
2588 copy_v3_v3(bp->prevpos, bp->pos);
2589 copy_v3_v3(bp->pos, bp->origT);
2590 }
2591 }
2592 }
2593
apply_spring_memory(Object * ob)2594 static void apply_spring_memory(Object *ob)
2595 {
2596 SoftBody *sb = ob->soft;
2597 BodySpring *bs;
2598 BodyPoint *bp1, *bp2;
2599 int a;
2600 float b, l, r;
2601
2602 if (sb && sb->totspring) {
2603 b = sb->plastic;
2604 for (a = 0; a < sb->totspring; a++) {
2605 bs = &sb->bspring[a];
2606 bp1 = &sb->bpoint[bs->v1];
2607 bp2 = &sb->bpoint[bs->v2];
2608 l = len_v3v3(bp1->pos, bp2->pos);
2609 r = bs->len / l;
2610 if ((r > 1.05f) || (r < 0.95f)) {
2611 bs->len = ((100.0f - b) * bs->len + b * l) / 100.0f;
2612 }
2613 }
2614 }
2615 }
2616
2617 /* expects full initialized softbody */
interpolate_exciter(Object * ob,int timescale,int time)2618 static void interpolate_exciter(Object *ob, int timescale, int time)
2619 {
2620 SoftBody *sb = ob->soft;
2621 BodyPoint *bp;
2622 float f;
2623 int a;
2624
2625 f = (float)time / (float)timescale;
2626
2627 for (a = sb->totpoint, bp = sb->bpoint; a > 0; a--, bp++) {
2628 bp->origT[0] = bp->origS[0] + f * (bp->origE[0] - bp->origS[0]);
2629 bp->origT[1] = bp->origS[1] + f * (bp->origE[1] - bp->origS[1]);
2630 bp->origT[2] = bp->origS[2] + f * (bp->origE[2] - bp->origS[2]);
2631 if (_final_goal(ob, bp) >= SOFTGOALSNAP) {
2632 bp->vec[0] = bp->origE[0] - bp->origS[0];
2633 bp->vec[1] = bp->origE[1] - bp->origS[1];
2634 bp->vec[2] = bp->origE[2] - bp->origS[2];
2635 }
2636 }
2637 }
2638
2639 /* ************ convertors ********** */
2640
2641 /* for each object type we need;
2642 * - xxxx_to_softbody(Object *ob) : a full (new) copy, creates SB geometry
2643 */
2644
2645 /* Resetting a Mesh SB object's springs */
2646 /* Spring length are caculted from'raw' mesh vertices that are NOT altered by modifier stack. */
springs_from_mesh(Object * ob)2647 static void springs_from_mesh(Object *ob)
2648 {
2649 SoftBody *sb;
2650 Mesh *me = ob->data;
2651 BodyPoint *bp;
2652 int a;
2653 float scale = 1.0f;
2654
2655 sb = ob->soft;
2656 if (me && sb) {
2657 /* using bp->origS as a container for spring calculations here
2658 * will be overwritten sbObjectStep() to receive
2659 * actual modifier stack positions
2660 */
2661 if (me->totvert) {
2662 bp = ob->soft->bpoint;
2663 for (a = 0; a < me->totvert; a++, bp++) {
2664 copy_v3_v3(bp->origS, me->mvert[a].co);
2665 mul_m4_v3(ob->obmat, bp->origS);
2666 }
2667 }
2668 /* recalculate spring length for meshes here */
2669 /* public version shrink to fit */
2670 if (sb->springpreload != 0) {
2671 scale = sb->springpreload / 100.0f;
2672 }
2673 for (a = 0; a < sb->totspring; a++) {
2674 BodySpring *bs = &sb->bspring[a];
2675 bs->len = scale * len_v3v3(sb->bpoint[bs->v1].origS, sb->bpoint[bs->v2].origS);
2676 }
2677 }
2678 }
2679
2680 /* makes totally fresh start situation */
mesh_to_softbody(Scene * scene,Object * ob)2681 static void mesh_to_softbody(Scene *scene, Object *ob)
2682 {
2683 SoftBody *sb;
2684 Mesh *me = ob->data;
2685 MEdge *medge = me->medge;
2686 BodyPoint *bp;
2687 BodySpring *bs;
2688 int a, totedge;
2689 int defgroup_index, defgroup_index_mass, defgroup_index_spring;
2690
2691 if (ob->softflag & OB_SB_EDGES) {
2692 totedge = me->totedge;
2693 }
2694 else {
2695 totedge = 0;
2696 }
2697
2698 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
2699 renew_softbody(scene, ob, me->totvert, totedge);
2700
2701 /* we always make body points */
2702 sb = ob->soft;
2703 bp = sb->bpoint;
2704
2705 defgroup_index = me->dvert ? (sb->vertgroup - 1) : -1;
2706 defgroup_index_mass = me->dvert ? BKE_object_defgroup_name_index(ob, sb->namedVG_Mass) : -1;
2707 defgroup_index_spring = me->dvert ? BKE_object_defgroup_name_index(ob, sb->namedVG_Spring_K) :
2708 -1;
2709
2710 for (a = 0; a < me->totvert; a++, bp++) {
2711 /* get scalar values needed *per vertex* from vertex group functions,
2712 * so we can *paint* them nicely ..
2713 * they are normalized [0.0..1.0] so may be we need amplitude for scale
2714 * which can be done by caller but still .. i'd like it to go this way
2715 */
2716
2717 if (ob->softflag & OB_SB_GOAL) {
2718 BLI_assert(bp->goal == sb->defgoal);
2719 }
2720 if ((ob->softflag & OB_SB_GOAL) && (defgroup_index != -1)) {
2721 bp->goal *= BKE_defvert_find_weight(&me->dvert[a], defgroup_index);
2722 }
2723
2724 /* to proof the concept
2725 * this enables per vertex *mass painting*
2726 */
2727
2728 if (defgroup_index_mass != -1) {
2729 bp->mass *= BKE_defvert_find_weight(&me->dvert[a], defgroup_index_mass);
2730 }
2731
2732 if (defgroup_index_spring != -1) {
2733 bp->springweight *= BKE_defvert_find_weight(&me->dvert[a], defgroup_index_spring);
2734 }
2735 }
2736
2737 /* but we only optionally add body edge springs */
2738 if (ob->softflag & OB_SB_EDGES) {
2739 if (medge) {
2740 bs = sb->bspring;
2741 for (a = me->totedge; a > 0; a--, medge++, bs++) {
2742 bs->v1 = medge->v1;
2743 bs->v2 = medge->v2;
2744 bs->springtype = SB_EDGE;
2745 }
2746
2747 /* insert *diagonal* springs in quads if desired */
2748 if (ob->softflag & OB_SB_QUADS) {
2749 add_mesh_quad_diag_springs(ob);
2750 }
2751
2752 build_bps_springlist(ob); /* scan for springs attached to bodypoints ONCE */
2753 /* insert *other second order* springs if desired */
2754 if (sb->secondspring > 0.0000001f) {
2755 /* exploits the first run of build_bps_springlist(ob); */
2756 add_2nd_order_springs(ob, sb->secondspring);
2757 /* yes we need to do it again. */
2758 build_bps_springlist(ob);
2759 }
2760 springs_from_mesh(ob); /* write the 'rest'-length of the springs */
2761 if (ob->softflag & OB_SB_SELF) {
2762 calculate_collision_balls(ob);
2763 }
2764 }
2765 }
2766 }
2767
mesh_faces_to_scratch(Object * ob)2768 static void mesh_faces_to_scratch(Object *ob)
2769 {
2770 SoftBody *sb = ob->soft;
2771 const Mesh *me = ob->data;
2772 MLoopTri *looptri, *lt;
2773 BodyFace *bodyface;
2774 int a;
2775 /* alloc and copy faces*/
2776
2777 sb->scratch->totface = poly_to_tri_count(me->totpoly, me->totloop);
2778 looptri = lt = MEM_mallocN(sizeof(*looptri) * sb->scratch->totface, __func__);
2779 BKE_mesh_recalc_looptri(me->mloop, me->mpoly, me->mvert, me->totloop, me->totpoly, looptri);
2780
2781 bodyface = sb->scratch->bodyface = MEM_mallocN(sizeof(BodyFace) * sb->scratch->totface,
2782 "SB_body_Faces");
2783
2784 for (a = 0; a < sb->scratch->totface; a++, lt++, bodyface++) {
2785 bodyface->v1 = me->mloop[lt->tri[0]].v;
2786 bodyface->v2 = me->mloop[lt->tri[1]].v;
2787 bodyface->v3 = me->mloop[lt->tri[2]].v;
2788 zero_v3(bodyface->ext_force);
2789 bodyface->ext_force[0] = bodyface->ext_force[1] = bodyface->ext_force[2] = 0.0f;
2790 bodyface->flag = 0;
2791 }
2792
2793 MEM_freeN(looptri);
2794 }
reference_to_scratch(Object * ob)2795 static void reference_to_scratch(Object *ob)
2796 {
2797 SoftBody *sb = ob->soft;
2798 ReferenceVert *rp;
2799 BodyPoint *bp;
2800 float accu_pos[3] = {0.f, 0.f, 0.f};
2801 float accu_mass = 0.f;
2802 int a;
2803
2804 sb->scratch->Ref.ivert = MEM_mallocN(sizeof(ReferenceVert) * sb->totpoint, "SB_Reference");
2805 bp = ob->soft->bpoint;
2806 rp = sb->scratch->Ref.ivert;
2807 for (a = 0; a < sb->totpoint; a++, rp++, bp++) {
2808 copy_v3_v3(rp->pos, bp->pos);
2809 add_v3_v3(accu_pos, bp->pos);
2810 accu_mass += _final_mass(ob, bp);
2811 }
2812 mul_v3_fl(accu_pos, 1.0f / accu_mass);
2813 copy_v3_v3(sb->scratch->Ref.com, accu_pos);
2814 /* printf("reference_to_scratch\n"); */
2815 }
2816
2817 /*
2818 * helper function to get proper spring length
2819 * when object is rescaled
2820 */
globallen(float * v1,float * v2,Object * ob)2821 static float globallen(float *v1, float *v2, Object *ob)
2822 {
2823 float p1[3], p2[3];
2824 copy_v3_v3(p1, v1);
2825 mul_m4_v3(ob->obmat, p1);
2826 copy_v3_v3(p2, v2);
2827 mul_m4_v3(ob->obmat, p2);
2828 return len_v3v3(p1, p2);
2829 }
2830
makelatticesprings(Lattice * lt,BodySpring * bs,int dostiff,Object * ob)2831 static void makelatticesprings(Lattice *lt, BodySpring *bs, int dostiff, Object *ob)
2832 {
2833 BPoint *bp = lt->def, *bpu;
2834 int u, v, w, dv, dw, bpc = 0, bpuc;
2835
2836 dv = lt->pntsu;
2837 dw = dv * lt->pntsv;
2838
2839 for (w = 0; w < lt->pntsw; w++) {
2840
2841 for (v = 0; v < lt->pntsv; v++) {
2842
2843 for (u = 0, bpuc = 0, bpu = NULL; u < lt->pntsu; u++, bp++, bpc++) {
2844
2845 if (w) {
2846 bs->v1 = bpc;
2847 bs->v2 = bpc - dw;
2848 bs->springtype = SB_EDGE;
2849 bs->len = globallen((bp - dw)->vec, bp->vec, ob);
2850 bs++;
2851 }
2852 if (v) {
2853 bs->v1 = bpc;
2854 bs->v2 = bpc - dv;
2855 bs->springtype = SB_EDGE;
2856 bs->len = globallen((bp - dv)->vec, bp->vec, ob);
2857 bs++;
2858 }
2859 if (u) {
2860 bs->v1 = bpuc;
2861 bs->v2 = bpc;
2862 bs->springtype = SB_EDGE;
2863 bs->len = globallen((bpu)->vec, bp->vec, ob);
2864 bs++;
2865 }
2866
2867 if (dostiff) {
2868
2869 if (w) {
2870 if (v && u) {
2871 bs->v1 = bpc;
2872 bs->v2 = bpc - dw - dv - 1;
2873 bs->springtype = SB_BEND;
2874 bs->len = globallen((bp - dw - dv - 1)->vec, bp->vec, ob);
2875 bs++;
2876 }
2877 if ((v < lt->pntsv - 1) && (u != 0)) {
2878 bs->v1 = bpc;
2879 bs->v2 = bpc - dw + dv - 1;
2880 bs->springtype = SB_BEND;
2881 bs->len = globallen((bp - dw + dv - 1)->vec, bp->vec, ob);
2882 bs++;
2883 }
2884 }
2885
2886 if (w < lt->pntsw - 1) {
2887 if (v && u) {
2888 bs->v1 = bpc;
2889 bs->v2 = bpc + dw - dv - 1;
2890 bs->springtype = SB_BEND;
2891 bs->len = globallen((bp + dw - dv - 1)->vec, bp->vec, ob);
2892 bs++;
2893 }
2894 if ((v < lt->pntsv - 1) && (u != 0)) {
2895 bs->v1 = bpc;
2896 bs->v2 = bpc + dw + dv - 1;
2897 bs->springtype = SB_BEND;
2898 bs->len = globallen((bp + dw + dv - 1)->vec, bp->vec, ob);
2899 bs++;
2900 }
2901 }
2902 }
2903 bpu = bp;
2904 bpuc = bpc;
2905 }
2906 }
2907 }
2908 }
2909
2910 /* makes totally fresh start situation */
lattice_to_softbody(Scene * scene,Object * ob)2911 static void lattice_to_softbody(Scene *scene, Object *ob)
2912 {
2913 Lattice *lt = ob->data;
2914 SoftBody *sb;
2915 int totvert, totspring = 0, a;
2916 BodyPoint *bp;
2917 BPoint *bpnt = lt->def;
2918 int defgroup_index, defgroup_index_mass, defgroup_index_spring;
2919
2920 totvert = lt->pntsu * lt->pntsv * lt->pntsw;
2921
2922 if (ob->softflag & OB_SB_EDGES) {
2923 totspring = ((lt->pntsu - 1) * lt->pntsv + (lt->pntsv - 1) * lt->pntsu) * lt->pntsw +
2924 lt->pntsu * lt->pntsv * (lt->pntsw - 1);
2925 if (ob->softflag & OB_SB_QUADS) {
2926 totspring += 4 * (lt->pntsu - 1) * (lt->pntsv - 1) * (lt->pntsw - 1);
2927 }
2928 }
2929
2930 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
2931 renew_softbody(scene, ob, totvert, totspring);
2932 sb = ob->soft; /* can be created in renew_softbody() */
2933 bp = sb->bpoint;
2934
2935 defgroup_index = lt->dvert ? (sb->vertgroup - 1) : -1;
2936 defgroup_index_mass = lt->dvert ? BKE_object_defgroup_name_index(ob, sb->namedVG_Mass) : -1;
2937 defgroup_index_spring = lt->dvert ? BKE_object_defgroup_name_index(ob, sb->namedVG_Spring_K) :
2938 -1;
2939
2940 /* same code used as for mesh vertices */
2941 for (a = 0; a < totvert; a++, bp++, bpnt++) {
2942
2943 if (ob->softflag & OB_SB_GOAL) {
2944 BLI_assert(bp->goal == sb->defgoal);
2945 }
2946
2947 if ((ob->softflag & OB_SB_GOAL) && (defgroup_index != -1)) {
2948 bp->goal *= BKE_defvert_find_weight(<->dvert[a], defgroup_index);
2949 }
2950 else {
2951 bp->goal *= bpnt->weight;
2952 }
2953
2954 if (defgroup_index_mass != -1) {
2955 bp->mass *= BKE_defvert_find_weight(<->dvert[a], defgroup_index_mass);
2956 }
2957
2958 if (defgroup_index_spring != -1) {
2959 bp->springweight *= BKE_defvert_find_weight(<->dvert[a], defgroup_index_spring);
2960 }
2961 }
2962
2963 /* create some helper edges to enable SB lattice to be useful at all */
2964 if (ob->softflag & OB_SB_EDGES) {
2965 makelatticesprings(lt, ob->soft->bspring, ob->softflag & OB_SB_QUADS, ob);
2966 build_bps_springlist(ob); /* link bps to springs */
2967 if (ob->softflag & OB_SB_SELF) {
2968 calculate_collision_balls(ob);
2969 }
2970 }
2971 }
2972
2973 /* makes totally fresh start situation */
curve_surf_to_softbody(Scene * scene,Object * ob)2974 static void curve_surf_to_softbody(Scene *scene, Object *ob)
2975 {
2976 Curve *cu = ob->data;
2977 SoftBody *sb;
2978 BodyPoint *bp;
2979 BodySpring *bs;
2980 Nurb *nu;
2981 BezTriple *bezt;
2982 BPoint *bpnt;
2983 int a, curindex = 0;
2984 int totvert, totspring = 0, setgoal = 0;
2985
2986 totvert = BKE_nurbList_verts_count(&cu->nurb);
2987
2988 if (ob->softflag & OB_SB_EDGES) {
2989 if (ob->type == OB_CURVE) {
2990 totspring = totvert - BLI_listbase_count(&cu->nurb);
2991 }
2992 }
2993
2994 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */
2995 renew_softbody(scene, ob, totvert, totspring);
2996 sb = ob->soft; /* can be created in renew_softbody() */
2997
2998 /* set vars now */
2999 bp = sb->bpoint;
3000 bs = sb->bspring;
3001
3002 /* weights from bpoints, same code used as for mesh vertices */
3003 /* if ((ob->softflag & OB_SB_GOAL) && sb->vertgroup) 2.4x hack*/
3004 /* new! take the weights from curve vertex anyhow */
3005 if (ob->softflag & OB_SB_GOAL) {
3006 setgoal = 1;
3007 }
3008
3009 for (nu = cu->nurb.first; nu; nu = nu->next) {
3010 if (nu->bezt) {
3011 /* Bezier case; this is nicly said naive; who ever wrote this part,
3012 * it was not me (JOW) :).
3013 *
3014 * a: never ever make tangent handles (sub) and or (ob)ject to collision.
3015 * b: rather calculate them using some C2
3016 * (C2= continuous in second derivative -> no jump in bending ) condition.
3017 *
3018 * Not too hard to do, but needs some more code to care for;
3019 * some one may want look at it (JOW 2010/06/12). */
3020 for (bezt = nu->bezt, a = 0; a < nu->pntsu; a++, bezt++, bp += 3, curindex += 3) {
3021 if (setgoal) {
3022 bp->goal *= bezt->weight;
3023
3024 /* all three triples */
3025 (bp + 1)->goal = bp->goal;
3026 (bp + 2)->goal = bp->goal;
3027 /*do not collide handles */
3028 (bp + 1)->loc_flag |= SBF_OUTOFCOLLISION;
3029 (bp + 2)->loc_flag |= SBF_OUTOFCOLLISION;
3030 }
3031
3032 if (totspring) {
3033 if (a > 0) {
3034 bs->v1 = curindex - 3;
3035 bs->v2 = curindex;
3036 bs->springtype = SB_HANDLE;
3037 bs->len = globallen((bezt - 1)->vec[0], bezt->vec[0], ob);
3038 bs++;
3039 }
3040 bs->v1 = curindex;
3041 bs->v2 = curindex + 1;
3042 bs->springtype = SB_HANDLE;
3043 bs->len = globallen(bezt->vec[0], bezt->vec[1], ob);
3044 bs++;
3045
3046 bs->v1 = curindex + 1;
3047 bs->v2 = curindex + 2;
3048 bs->springtype = SB_HANDLE;
3049 bs->len = globallen(bezt->vec[1], bezt->vec[2], ob);
3050 bs++;
3051 }
3052 }
3053 }
3054 else {
3055 for (bpnt = nu->bp, a = 0; a < nu->pntsu * nu->pntsv; a++, bpnt++, bp++, curindex++) {
3056 if (setgoal) {
3057 bp->goal *= bpnt->weight;
3058 }
3059 if (totspring && a > 0) {
3060 bs->v1 = curindex - 1;
3061 bs->v2 = curindex;
3062 bs->springtype = SB_EDGE;
3063 bs->len = globallen((bpnt - 1)->vec, bpnt->vec, ob);
3064 bs++;
3065 }
3066 }
3067 }
3068 }
3069
3070 if (totspring) {
3071 build_bps_springlist(ob); /* link bps to springs */
3072 if (ob->softflag & OB_SB_SELF) {
3073 calculate_collision_balls(ob);
3074 }
3075 }
3076 }
3077
3078 /* copies softbody result back in object */
softbody_to_object(Object * ob,float (* vertexCos)[3],int numVerts,int local)3079 static void softbody_to_object(Object *ob, float (*vertexCos)[3], int numVerts, int local)
3080 {
3081 SoftBody *sb = ob->soft;
3082 if (sb) {
3083 BodyPoint *bp = sb->bpoint;
3084 int a;
3085 if (sb->solverflags & SBSO_ESTIMATEIPO) {
3086 SB_estimate_transform(ob, sb->lcom, sb->lrot, sb->lscale);
3087 }
3088 /* inverse matrix is not uptodate... */
3089 invert_m4_m4(ob->imat, ob->obmat);
3090
3091 for (a = 0; a < numVerts; a++, bp++) {
3092 copy_v3_v3(vertexCos[a], bp->pos);
3093 if (local == 0) {
3094 mul_m4_v3(ob->imat, vertexCos[a]); /* softbody is in global coords, baked optionally not */
3095 }
3096 }
3097 }
3098 }
3099
3100 /* +++ ************ maintaining scratch *************** */
sb_new_scratch(SoftBody * sb)3101 static void sb_new_scratch(SoftBody *sb)
3102 {
3103 if (!sb) {
3104 return;
3105 }
3106 sb->scratch = MEM_callocN(sizeof(SBScratch), "SBScratch");
3107 sb->scratch->colliderhash = BLI_ghash_ptr_new("sb_new_scratch gh");
3108 sb->scratch->bodyface = NULL;
3109 sb->scratch->totface = 0;
3110 sb->scratch->aabbmax[0] = sb->scratch->aabbmax[1] = sb->scratch->aabbmax[2] = 1.0e30f;
3111 sb->scratch->aabbmin[0] = sb->scratch->aabbmin[1] = sb->scratch->aabbmin[2] = -1.0e30f;
3112 sb->scratch->Ref.ivert = NULL;
3113 }
3114 /* --- ************ maintaining scratch *************** */
3115
3116 /* ************ Object level, exported functions *************** */
3117
3118 /* allocates and initializes general main data */
sbNew(Scene * scene)3119 SoftBody *sbNew(Scene *scene)
3120 {
3121 SoftBody *sb;
3122
3123 sb = MEM_callocN(sizeof(SoftBody), "softbody");
3124
3125 sb->mediafrict = 0.5f;
3126 sb->nodemass = 1.0f;
3127 sb->grav = 9.8f;
3128 sb->physics_speed = 1.0f;
3129 sb->rklimit = 0.1f;
3130
3131 sb->goalspring = 0.5f;
3132 sb->goalfrict = 0.0f;
3133 sb->mingoal = 0.0f;
3134 sb->maxgoal = 1.0f;
3135 sb->defgoal = 0.7f;
3136
3137 sb->inspring = 0.5f;
3138 sb->infrict = 0.5f;
3139 /*todo backward file compat should copy inspring to inpush while reading old files*/
3140 sb->inpush = 0.5f;
3141
3142 sb->interval = 10;
3143 if (scene != NULL) {
3144 sb->sfra = scene->r.sfra;
3145 sb->efra = scene->r.efra;
3146 }
3147
3148 sb->colball = 0.49f;
3149 sb->balldamp = 0.50f;
3150 sb->ballstiff = 1.0f;
3151 sb->sbc_mode = 1;
3152
3153 sb->minloops = 10;
3154 sb->maxloops = 300;
3155
3156 sb->choke = 3;
3157 sb_new_scratch(sb);
3158 /*todo backward file compat should set sb->shearstiff = 1.0f while reading old files*/
3159 sb->shearstiff = 1.0f;
3160 sb->solverflags |= SBSO_OLDERR;
3161
3162 sb->shared = MEM_callocN(sizeof(*sb->shared), "SoftBody_Shared");
3163 sb->shared->pointcache = BKE_ptcache_add(&sb->shared->ptcaches);
3164
3165 if (!sb->effector_weights) {
3166 sb->effector_weights = BKE_effector_add_weights(NULL);
3167 }
3168
3169 sb->last_frame = MINFRAME - 1;
3170
3171 return sb;
3172 }
3173
3174 /* frees all */
sbFree(Object * ob)3175 void sbFree(Object *ob)
3176 {
3177 SoftBody *sb = ob->soft;
3178 if (sb == NULL) {
3179 return;
3180 }
3181
3182 free_softbody_intern(sb);
3183
3184 if ((ob->id.tag & LIB_TAG_NO_MAIN) == 0) {
3185 /* Only free shared data on non-CoW copies */
3186 BKE_ptcache_free_list(&sb->shared->ptcaches);
3187 sb->shared->pointcache = NULL;
3188 MEM_freeN(sb->shared);
3189 }
3190 if (sb->effector_weights) {
3191 MEM_freeN(sb->effector_weights);
3192 }
3193 MEM_freeN(sb);
3194
3195 ob->soft = NULL;
3196 }
3197
sbFreeSimulation(SoftBody * sb)3198 void sbFreeSimulation(SoftBody *sb)
3199 {
3200 free_softbody_intern(sb);
3201 }
3202
3203 /* makes totally fresh start situation */
sbObjectToSoftbody(Object * ob)3204 void sbObjectToSoftbody(Object *ob)
3205 {
3206 // ob->softflag |= OB_SB_REDO;
3207
3208 free_softbody_intern(ob->soft);
3209 }
3210
object_has_edges(Object * ob)3211 static bool object_has_edges(Object *ob)
3212 {
3213 if (ob->type == OB_MESH) {
3214 return ((Mesh *)ob->data)->totedge;
3215 }
3216 if (ob->type == OB_LATTICE) {
3217 return true;
3218 }
3219
3220 return false;
3221 }
3222
3223 /* SB global visible functions */
sbSetInterruptCallBack(int (* f)(void))3224 void sbSetInterruptCallBack(int (*f)(void))
3225 {
3226 SB_localInterruptCallBack = f;
3227 }
3228
softbody_update_positions(Object * ob,SoftBody * sb,float (* vertexCos)[3],int numVerts)3229 static void softbody_update_positions(Object *ob,
3230 SoftBody *sb,
3231 float (*vertexCos)[3],
3232 int numVerts)
3233 {
3234 BodyPoint *bp;
3235 int a;
3236
3237 if (!sb || !sb->bpoint) {
3238 return;
3239 }
3240
3241 for (a = 0, bp = sb->bpoint; a < numVerts; a++, bp++) {
3242 /* store where goals are now */
3243 copy_v3_v3(bp->origS, bp->origE);
3244 /* copy the position of the goals at desired end time */
3245 copy_v3_v3(bp->origE, vertexCos[a]);
3246 /* vertexCos came from local world, go global */
3247 mul_m4_v3(ob->obmat, bp->origE);
3248 /* just to be save give bp->origT a defined value
3249 * will be calculated in interpolate_exciter() */
3250 copy_v3_v3(bp->origT, bp->origE);
3251 }
3252 }
3253
3254 /* void SB_estimate_transform */
3255 /* input Object *ob out (says any object that can do SB like mesh, lattice, curve )
3256 * output float lloc[3], float lrot[3][3], float lscale[3][3]
3257 * that is:
3258 * a precise position vector denoting the motion of the center of mass
3259 * give a rotation/scale matrix using averaging method, that's why estimate and not calculate
3260 * see: this is kind of reverse engineering: having to states of a point cloud and recover what
3261 * happened our advantage here we know the identity of the vertex there are others methods giving
3262 * other results. lloc, lrot, lscale are allowed to be NULL, just in case you don't need it.
3263 * should be pretty useful for pythoneers :)
3264 * not! velocity .. 2nd order stuff
3265 * vcloud_estimate_transform_v3 see
3266 */
3267
SB_estimate_transform(Object * ob,float lloc[3],float lrot[3][3],float lscale[3][3])3268 void SB_estimate_transform(Object *ob, float lloc[3], float lrot[3][3], float lscale[3][3])
3269 {
3270 BodyPoint *bp;
3271 ReferenceVert *rp;
3272 SoftBody *sb = NULL;
3273 float(*opos)[3];
3274 float(*rpos)[3];
3275 float com[3], rcom[3];
3276 int a;
3277
3278 if (!ob || !ob->soft) {
3279 return; /* why did we get here ? */
3280 }
3281 sb = ob->soft;
3282 if (!sb || !sb->bpoint) {
3283 return;
3284 }
3285 opos = MEM_callocN(sizeof(float[3]) * sb->totpoint, "SB_OPOS");
3286 rpos = MEM_callocN(sizeof(float[3]) * sb->totpoint, "SB_RPOS");
3287 /* might filter vertex selection with a vertex group */
3288 for (a = 0, bp = sb->bpoint, rp = sb->scratch->Ref.ivert; a < sb->totpoint; a++, bp++, rp++) {
3289 copy_v3_v3(rpos[a], rp->pos);
3290 copy_v3_v3(opos[a], bp->pos);
3291 }
3292
3293 vcloud_estimate_transform_v3(sb->totpoint, opos, NULL, rpos, NULL, com, rcom, lrot, lscale);
3294 // sub_v3_v3(com, rcom);
3295 if (lloc) {
3296 copy_v3_v3(lloc, com);
3297 }
3298 copy_v3_v3(sb->lcom, com);
3299 if (lscale) {
3300 copy_m3_m3(sb->lscale, lscale);
3301 }
3302 if (lrot) {
3303 copy_m3_m3(sb->lrot, lrot);
3304 }
3305
3306 MEM_freeN(opos);
3307 MEM_freeN(rpos);
3308 }
3309
softbody_reset(Object * ob,SoftBody * sb,float (* vertexCos)[3],int numVerts)3310 static void softbody_reset(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts)
3311 {
3312 BodyPoint *bp;
3313 int a;
3314
3315 for (a = 0, bp = sb->bpoint; a < numVerts; a++, bp++) {
3316 copy_v3_v3(bp->pos, vertexCos[a]);
3317 mul_m4_v3(ob->obmat, bp->pos); /* yep, sofbody is global coords*/
3318 copy_v3_v3(bp->origS, bp->pos);
3319 copy_v3_v3(bp->origE, bp->pos);
3320 copy_v3_v3(bp->origT, bp->pos);
3321 bp->vec[0] = bp->vec[1] = bp->vec[2] = 0.0f;
3322
3323 /* the bp->prev*'s are for rolling back from a canceled try to propagate in time
3324 * adaptive step size algo in a nutshell:
3325 * 1. set scheduled time step to new dtime
3326 * 2. try to advance the scheduled time step, being optimistic execute it
3327 * 3. check for success
3328 * 3.a we 're fine continue, may be we can increase scheduled time again ?? if so, do so!
3329 * 3.b we did exceed error limit --> roll back, shorten the scheduled time and try again at 2.
3330 * 4. check if we did reach dtime
3331 * 4.a nope we need to do some more at 2.
3332 * 4.b yup we're done
3333 */
3334
3335 copy_v3_v3(bp->prevpos, bp->pos);
3336 copy_v3_v3(bp->prevvec, bp->vec);
3337 copy_v3_v3(bp->prevdx, bp->vec);
3338 copy_v3_v3(bp->prevdv, bp->vec);
3339 }
3340
3341 /* make a nice clean scratch struct */
3342 free_scratch(sb); /* clear if any */
3343 sb_new_scratch(sb); /* make a new */
3344 sb->scratch->needstobuildcollider = 1;
3345 zero_v3(sb->lcom);
3346 unit_m3(sb->lrot);
3347 unit_m3(sb->lscale);
3348
3349 /* copy some info to scratch */
3350 /* we only need that if we want to reconstruct IPO */
3351 if (1) {
3352 reference_to_scratch(ob);
3353 SB_estimate_transform(ob, NULL, NULL, NULL);
3354 SB_estimate_transform(ob, NULL, NULL, NULL);
3355 }
3356 switch (ob->type) {
3357 case OB_MESH:
3358 if (ob->softflag & OB_SB_FACECOLL) {
3359 mesh_faces_to_scratch(ob);
3360 }
3361 break;
3362 case OB_LATTICE:
3363 break;
3364 case OB_CURVE:
3365 case OB_SURF:
3366 break;
3367 default:
3368 break;
3369 }
3370 }
3371
softbody_step(struct Depsgraph * depsgraph,Scene * scene,Object * ob,SoftBody * sb,float dtime)3372 static void softbody_step(
3373 struct Depsgraph *depsgraph, Scene *scene, Object *ob, SoftBody *sb, float dtime)
3374 {
3375 /* the simulator */
3376 float forcetime;
3377 double sct, sst;
3378
3379 sst = PIL_check_seconds_timer();
3380 /* Integration back in time is possible in theory, but pretty useless here.
3381 * So we refuse to do so. Since we do not know anything about 'outside' changes
3382 * especially colliders we refuse to go more than 10 frames.
3383 */
3384 if (dtime < 0 || dtime > 10.5f) {
3385 return;
3386 }
3387
3388 ccd_update_deflector_hash(depsgraph, sb->collision_group, ob, sb->scratch->colliderhash);
3389
3390 if (sb->scratch->needstobuildcollider) {
3391 ccd_build_deflector_hash(depsgraph, sb->collision_group, ob, sb->scratch->colliderhash);
3392 sb->scratch->needstobuildcollider = 0;
3393 }
3394
3395 if (sb->solver_ID < 2) {
3396 /* special case of 2nd order Runge-Kutta type AKA Heun */
3397 int mid_flags = 0;
3398 float err = 0;
3399 /* Set defaults guess we shall do one frame */
3400 float forcetimemax = 1.0f;
3401 /* Set defaults guess 1/100 is tight enough */
3402 float forcetimemin = 0.01f;
3403 /* How far did we get without violating error condition. */
3404 float timedone = 0.0;
3405 /* Loops = counter for emergency brake we don't want to lock up the system if physics fail. */
3406 int loops = 0;
3407
3408 SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */
3409 /* adjust loop limits */
3410 if (sb->minloops > 0) {
3411 forcetimemax = dtime / sb->minloops;
3412 }
3413 if (sb->maxloops > 0) {
3414 forcetimemin = dtime / sb->maxloops;
3415 }
3416
3417 if (sb->solver_ID > 0) {
3418 mid_flags |= MID_PRESERVE;
3419 }
3420
3421 forcetime = forcetimemax; /* hope for integrating in one step */
3422 while ((fabsf(timedone) < fabsf(dtime)) && (loops < 2000)) {
3423 /* set goals in time */
3424 interpolate_exciter(ob, 200, (int)(200.0f * (timedone / dtime)));
3425
3426 sb->scratch->flag &= ~SBF_DOFUZZY;
3427 /* do predictive euler step */
3428 softbody_calc_forces(depsgraph, scene, ob, forcetime, timedone / dtime);
3429
3430 softbody_apply_forces(ob, forcetime, 1, NULL, mid_flags);
3431
3432 /* crop new slope values to do averaged slope step */
3433 softbody_calc_forces(depsgraph, scene, ob, forcetime, timedone / dtime);
3434
3435 softbody_apply_forces(ob, forcetime, 2, &err, mid_flags);
3436 softbody_apply_goalsnap(ob);
3437
3438 if (err > SoftHeunTol) { /* error needs to be scaled to some quantity */
3439
3440 if (forcetime > forcetimemin) {
3441 forcetime = max_ff(forcetime / 2.0f, forcetimemin);
3442 softbody_restore_prev_step(ob);
3443 // printf("down, ");
3444 }
3445 else {
3446 timedone += forcetime;
3447 }
3448 }
3449 else {
3450 float newtime = forcetime * 1.1f; /* hope for 1.1 times better conditions in next step */
3451
3452 if (sb->scratch->flag & SBF_DOFUZZY) {
3453 ///* stay with this stepsize unless err really small */
3454 // if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) {
3455 newtime = forcetime;
3456 //}
3457 }
3458 else {
3459 if (err > SoftHeunTol / 2.0f) { /* stay with this stepsize unless err really small */
3460 newtime = forcetime;
3461 }
3462 }
3463 timedone += forcetime;
3464 newtime = min_ff(forcetimemax, max_ff(newtime, forcetimemin));
3465 // if (newtime > forcetime) printf("up, ");
3466 if (forcetime > 0.0f) {
3467 forcetime = min_ff(dtime - timedone, newtime);
3468 }
3469 else {
3470 forcetime = max_ff(dtime - timedone, newtime);
3471 }
3472 }
3473 loops++;
3474 if (sb->solverflags & SBSO_MONITOR) {
3475 sct = PIL_check_seconds_timer();
3476 if (sct - sst > 0.5) {
3477 printf("%3.0f%% \r", 100.0f * timedone / dtime);
3478 }
3479 }
3480 /* ask for user break */
3481 if (SB_localInterruptCallBack && SB_localInterruptCallBack()) {
3482 break;
3483 }
3484 }
3485 /* move snapped to final position */
3486 interpolate_exciter(ob, 2, 2);
3487 softbody_apply_goalsnap(ob);
3488
3489 // if (G.debug & G_DEBUG) {
3490 if (sb->solverflags & SBSO_MONITOR) {
3491 if (loops > HEUNWARNLIMIT) { /* monitor high loop counts */
3492 printf("\r needed %d steps/frame", loops);
3493 }
3494 }
3495 }
3496 else if (sb->solver_ID == 2) {
3497 /* do semi "fake" implicit euler */
3498 /* removed */
3499 } /*SOLVER SELECT*/
3500 else if (sb->solver_ID == 4) {
3501 /* do semi "fake" implicit euler */
3502 } /*SOLVER SELECT*/
3503 else if (sb->solver_ID == 3) {
3504 /* do "stupid" semi "fake" implicit euler */
3505 /* removed */
3506
3507 } /*SOLVER SELECT*/
3508 else {
3509 CLOG_ERROR(&LOG, "softbody no valid solver ID!");
3510 } /*SOLVER SELECT*/
3511 if (sb->plastic) {
3512 apply_spring_memory(ob);
3513 }
3514
3515 if (sb->solverflags & SBSO_MONITOR) {
3516 sct = PIL_check_seconds_timer();
3517 if ((sct - sst > 0.5) || (G.debug & G_DEBUG)) {
3518 printf(" solver time %f sec %s\n", sct - sst, ob->id.name);
3519 }
3520 }
3521 }
3522
sbStoreLastFrame(struct Depsgraph * depsgraph,Object * object,float framenr)3523 static void sbStoreLastFrame(struct Depsgraph *depsgraph, Object *object, float framenr)
3524 {
3525 if (!DEG_is_active(depsgraph)) {
3526 return;
3527 }
3528 Object *object_orig = DEG_get_original_object(object);
3529 object->soft->last_frame = framenr;
3530 object_orig->soft->last_frame = framenr;
3531 }
3532
3533 /* simulates one step. framenr is in frames */
sbObjectStep(struct Depsgraph * depsgraph,Scene * scene,Object * ob,float cfra,float (* vertexCos)[3],int numVerts)3534 void sbObjectStep(struct Depsgraph *depsgraph,
3535 Scene *scene,
3536 Object *ob,
3537 float cfra,
3538 float (*vertexCos)[3],
3539 int numVerts)
3540 {
3541 SoftBody *sb = ob->soft;
3542 PointCache *cache;
3543 PTCacheID pid;
3544 float dtime, timescale;
3545 int framedelta, framenr, startframe, endframe;
3546 int cache_result;
3547 cache = sb->shared->pointcache;
3548
3549 framenr = (int)cfra;
3550 framedelta = framenr - cache->simframe;
3551
3552 BKE_ptcache_id_from_softbody(&pid, ob, sb);
3553 BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, ×cale);
3554
3555 /* check for changes in mesh, should only happen in case the mesh
3556 * structure changes during an animation */
3557 if (sb->bpoint && numVerts != sb->totpoint) {
3558 BKE_ptcache_invalidate(cache);
3559 return;
3560 }
3561
3562 /* clamp frame ranges */
3563 if (framenr < startframe) {
3564 BKE_ptcache_invalidate(cache);
3565 return;
3566 }
3567 if (framenr > endframe) {
3568 framenr = endframe;
3569 }
3570
3571 /* verify if we need to create the softbody data */
3572 if (sb->bpoint == NULL ||
3573 ((ob->softflag & OB_SB_EDGES) && !ob->soft->bspring && object_has_edges(ob))) {
3574
3575 switch (ob->type) {
3576 case OB_MESH:
3577 mesh_to_softbody(scene, ob);
3578 break;
3579 case OB_LATTICE:
3580 lattice_to_softbody(scene, ob);
3581 break;
3582 case OB_CURVE:
3583 case OB_SURF:
3584 curve_surf_to_softbody(scene, ob);
3585 break;
3586 default:
3587 renew_softbody(scene, ob, numVerts, 0);
3588 break;
3589 }
3590
3591 softbody_update_positions(ob, sb, vertexCos, numVerts);
3592 softbody_reset(ob, sb, vertexCos, numVerts);
3593 }
3594
3595 /* still no points? go away */
3596 if (sb->totpoint == 0) {
3597 return;
3598 }
3599 if (framenr == startframe) {
3600 BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
3601
3602 /* first frame, no simulation to do, just set the positions */
3603 softbody_update_positions(ob, sb, vertexCos, numVerts);
3604
3605 BKE_ptcache_validate(cache, framenr);
3606 cache->flag &= ~PTCACHE_REDO_NEEDED;
3607
3608 sbStoreLastFrame(depsgraph, ob, framenr);
3609
3610 return;
3611 }
3612
3613 /* try to read from cache */
3614 bool can_write_cache = DEG_is_active(depsgraph);
3615 bool can_simulate = (framenr == sb->last_frame + 1) && !(cache->flag & PTCACHE_BAKED) &&
3616 can_write_cache;
3617
3618 cache_result = BKE_ptcache_read(&pid, (float)framenr + scene->r.subframe, can_simulate);
3619
3620 if (cache_result == PTCACHE_READ_EXACT || cache_result == PTCACHE_READ_INTERPOLATED ||
3621 (!can_simulate && cache_result == PTCACHE_READ_OLD)) {
3622 softbody_to_object(ob, vertexCos, numVerts, sb->local);
3623
3624 BKE_ptcache_validate(cache, framenr);
3625
3626 if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED &&
3627 can_write_cache) {
3628 BKE_ptcache_write(&pid, framenr);
3629 }
3630
3631 sbStoreLastFrame(depsgraph, ob, framenr);
3632
3633 return;
3634 }
3635 if (cache_result == PTCACHE_READ_OLD) {
3636 /* pass */
3637 }
3638 else if (/*ob->id.lib || */
3639 /* "library linking & pointcaches" has to be solved properly at some point */
3640 (cache->flag & PTCACHE_BAKED)) {
3641 /* if baked and nothing in cache, do nothing */
3642 if (can_write_cache) {
3643 BKE_ptcache_invalidate(cache);
3644 }
3645 return;
3646 }
3647
3648 if (!can_simulate) {
3649 return;
3650 }
3651
3652 /* if on second frame, write cache for first frame */
3653 if (cache->simframe == startframe &&
3654 (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
3655 BKE_ptcache_write(&pid, startframe);
3656 }
3657
3658 softbody_update_positions(ob, sb, vertexCos, numVerts);
3659
3660 /* checking time: */
3661 dtime = framedelta * timescale;
3662
3663 /* do simulation */
3664 softbody_step(depsgraph, scene, ob, sb, dtime);
3665
3666 softbody_to_object(ob, vertexCos, numVerts, 0);
3667
3668 BKE_ptcache_validate(cache, framenr);
3669 BKE_ptcache_write(&pid, framenr);
3670
3671 sbStoreLastFrame(depsgraph, ob, framenr);
3672 }
3673