|
Blender
V2.59
|
00001 /* softbody.c 00002 * 00003 * $Id: softbody.c 38235 2011-07-08 13:22:58Z blendix $ 00004 * 00005 * ***** BEGIN GPL LICENSE BLOCK ***** 00006 * 00007 * This program is free software; you can redistribute it and/or 00008 * modify it under the terms of the GNU General Public License 00009 * as published by the Free Software Foundation; either version 2 00010 * of the License, or (at your option) any later version. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU General Public License 00018 * along with this program; if not, write to the Free Software Foundation, 00019 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00020 * 00021 * The Original Code is Copyright (C) Blender Foundation 00022 * All rights reserved. 00023 * 00024 * The Original Code is: all of this file. 00025 * 00026 * Contributor(s): none yet. 00027 * 00028 * ***** END GPL LICENSE BLOCK ***** 00029 */ 00030 00036 /* 00037 ****** 00038 variables on the UI for now 00039 00040 float mediafrict; friction to env 00041 float nodemass; softbody mass of *vertex* 00042 float grav; softbody amount of gravitaion to apply 00043 00044 float goalspring; softbody goal springs 00045 float goalfrict; softbody goal springs friction 00046 float mingoal; quick limits for goal 00047 float maxgoal; 00048 00049 float inspring; softbody inner springs 00050 float infrict; softbody inner springs friction 00051 00052 ***** 00053 */ 00054 00055 00056 #include <math.h> 00057 #include <stdlib.h> 00058 #include <string.h> 00059 00060 #include "MEM_guardedalloc.h" 00061 00062 /* types */ 00063 #include "DNA_object_types.h" 00064 #include "DNA_scene_types.h" 00065 #include "DNA_lattice_types.h" 00066 #include "DNA_curve_types.h" 00067 #include "DNA_mesh_types.h" 00068 #include "DNA_meshdata_types.h" 00069 00070 #include "BLI_math.h" 00071 #include "BLI_utildefines.h" 00072 #include "BLI_ghash.h" 00073 #include "BLI_threads.h" 00074 00075 #include "BKE_curve.h" 00076 #include "BKE_effect.h" 00077 #include "BKE_global.h" 00078 #include "BKE_modifier.h" 00079 #include "BKE_softbody.h" 00080 #include "BKE_DerivedMesh.h" 00081 #include "BKE_pointcache.h" 00082 #include "BKE_deform.h" 00083 //XXX #include "BIF_editdeform.h" 00084 //XXX #include "BIF_graphics.h" 00085 #include "PIL_time.h" 00086 // #include "ONL_opennl.h" remove linking to ONL for now 00087 00088 /* callbacks for errors and interrupts and some goo */ 00089 static int (*SB_localInterruptCallBack)(void) = NULL; 00090 00091 00092 /* ********** soft body engine ******* */ 00093 00094 typedef enum {SB_EDGE=1,SB_BEND=2,SB_STIFFQUAD=3,SB_HANDLE=4} type_spring; 00095 00096 typedef struct BodySpring { 00097 int v1, v2; 00098 float len,cf,load; 00099 float ext_force[3]; /* edges colliding and sailing */ 00100 type_spring springtype; 00101 short flag; 00102 } BodySpring; 00103 00104 typedef struct BodyFace { 00105 int v1, v2, v3 ,v4; 00106 float ext_force[3]; /* faces colliding */ 00107 short flag; 00108 } BodyFace; 00109 00110 typedef struct ReferenceVert { 00111 float pos[3]; /* position relative to com */ 00112 float mass; /* node mass */ 00113 } ReferenceVert; 00114 00115 typedef struct ReferenceState { 00116 float com[3]; /* center of mass*/ 00117 ReferenceVert *ivert; /* list of intial values */ 00118 }ReferenceState ; 00119 00120 00121 /*private scratch pad for caching and other data only needed when alive*/ 00122 typedef struct SBScratch { 00123 GHash *colliderhash; 00124 short needstobuildcollider; 00125 short flag; 00126 BodyFace *bodyface; 00127 int totface; 00128 float aabbmin[3],aabbmax[3]; 00129 ReferenceState Ref; 00130 }SBScratch; 00131 00132 typedef struct SB_thread_context { 00133 Scene *scene; 00134 Object *ob; 00135 float forcetime; 00136 float timenow; 00137 int ifirst; 00138 int ilast; 00139 ListBase *do_effector; 00140 int do_deflector; 00141 float fieldfactor; 00142 float windfactor; 00143 int nr; 00144 int tot; 00145 }SB_thread_context; 00146 00147 #define NLF_BUILD 1 00148 #define NLF_SOLVE 2 00149 00150 #define MID_PRESERVE 1 00151 00152 #define SOFTGOALSNAP 0.999f 00153 /* if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp 00154 removes *unnecessary* stiffnes from ODE system 00155 */ 00156 #define HEUNWARNLIMIT 1 /* 500 would be fine i think for detecting severe *stiff* stuff */ 00157 00158 00159 #define BSF_INTERSECT 1 /* edge intersects collider face */ 00160 00161 /* private definitions for bodypoint states */ 00162 #define SBF_DOFUZZY 1 /* Bodypoint do fuzzy */ 00163 #define SBF_OUTOFCOLLISION 2 /* Bodypoint does not collide */ 00164 00165 00166 #define BFF_INTERSECT 1 /* collider edge intrudes face */ 00167 #define BFF_CLOSEVERT 2 /* collider vertex repulses face */ 00168 00169 00170 static float SoftHeunTol = 1.0f; /* humm .. this should be calculated from sb parameters and sizes */ 00171 00172 /* local prototypes */ 00173 static void free_softbody_intern(SoftBody *sb); 00174 /* aye this belongs to arith.c */ 00175 static void Vec3PlusStVec(float *v, float s, float *v1); 00176 00177 /*+++ frame based timing +++*/ 00178 00179 /*physical unit of force is [kg * m / sec^2]*/ 00180 00181 static float sb_grav_force_scale(Object *UNUSED(ob)) 00182 /* since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm 00183 put it to a function here, so we can add user options later without touching simulation code 00184 */ 00185 { 00186 return (0.001f); 00187 } 00188 00189 static float sb_fric_force_scale(Object *UNUSED(ob)) 00190 /* rescaling unit of drag [1 / sec] to somehow reasonable 00191 put it to a function here, so we can add user options later without touching simulation code 00192 */ 00193 { 00194 return (0.01f); 00195 } 00196 00197 static float sb_time_scale(Object *ob) 00198 /* defining the frames to *real* time relation */ 00199 { 00200 SoftBody *sb= ob->soft; /* is supposed to be there */ 00201 if (sb){ 00202 return(sb->physics_speed); 00203 /*hrms .. this could be IPO as well :) 00204 estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)] 00205 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions] has period 65 frames 00206 theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM) 00207 */ 00208 } 00209 return (1.0f); 00210 /* 00211 this would be frames/sec independant timing assuming 25 fps is default 00212 but does not work very well with NLA 00213 return (25.0f/scene->r.frs_sec) 00214 */ 00215 } 00216 /*--- frame based timing ---*/ 00217 00218 /* helper functions for everything is animatable jow_go_for2_5 +++++++*/ 00219 /* introducing them here, because i know: steps in properties ( at frame timing ) 00220 will cause unwanted responses of the softbody system (which does inter frame calculations ) 00221 so first 'cure' would be: interpolate linear in time .. 00222 Q: why do i write this? 00223 A: because it happend once, that some eger coder 'streamlined' code to fail. 00224 We DO linear interpolation for goals .. and i think we should do on animated properties as well 00225 */ 00226 00227 /* animate sb->maxgoal,sb->mingoal */ 00228 static float _final_goal(Object *ob,BodyPoint *bp)/*jow_go_for2_5 */ 00229 { 00230 float f = -1999.99f; 00231 if (ob){ 00232 SoftBody *sb= ob->soft; /* is supposed to be there */ 00233 if(!(ob->softflag & OB_SB_GOAL)) return (0.0f); 00234 if (sb&&bp){ 00235 if (bp->goal < 0.0f) return (0.0f); 00236 f = sb->mingoal + bp->goal*ABS(sb->maxgoal - sb->mingoal); 00237 f = pow(f, 4.0f); 00238 return (f); 00239 } 00240 } 00241 printf("_final_goal failed! sb or bp ==NULL \n" ); 00242 return f; /*using crude but spot able values some times helps debuggin */ 00243 } 00244 00245 static float _final_mass(Object *ob,BodyPoint *bp) 00246 { 00247 if (ob){ 00248 SoftBody *sb= ob->soft; /* is supposed to be there */ 00249 if (sb&&bp){ 00250 return(bp->mass*sb->nodemass); 00251 } 00252 } 00253 printf("_final_mass failed! sb or bp ==NULL \n" ); 00254 return 1.0f; 00255 } 00256 /* helper functions for everything is animateble jow_go_for2_5 ------*/ 00257 00258 /*+++ collider caching and dicing +++*/ 00259 00260 /******************** 00261 for each target object/face the axis aligned bounding box (AABB) is stored 00262 faces parallel to global axes 00263 so only simple "value" in [min,max] ckecks are used 00264 float operations still 00265 */ 00266 00267 /* just an ID here to reduce the prob for killing objects 00268 ** ob->sumohandle points to we should not kill :) 00269 */ 00270 static const int CCD_SAVETY = 190561; 00271 00272 typedef struct ccdf_minmax{ 00273 float minx,miny,minz,maxx,maxy,maxz; 00274 }ccdf_minmax; 00275 00276 00277 00278 typedef struct ccd_Mesh { 00279 int totvert, totface; 00280 MVert *mvert; 00281 MVert *mprevvert; 00282 MFace *mface; 00283 int savety; 00284 ccdf_minmax *mima; 00285 /* Axis Aligned Bounding Box AABB */ 00286 float bbmin[3]; 00287 float bbmax[3]; 00288 }ccd_Mesh; 00289 00290 00291 00292 00293 static ccd_Mesh *ccd_mesh_make(Object *ob) 00294 { 00295 CollisionModifierData *cmd; 00296 ccd_Mesh *pccd_M = NULL; 00297 ccdf_minmax *mima =NULL; 00298 MFace *mface=NULL; 00299 float v[3],hull; 00300 int i; 00301 00302 cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision); 00303 00304 /* first some paranoia checks */ 00305 if (!cmd) return NULL; 00306 if (!cmd->numverts || !cmd->numfaces) return NULL; 00307 00308 pccd_M = MEM_mallocN(sizeof(ccd_Mesh),"ccd_Mesh"); 00309 pccd_M->totvert = cmd->numverts; 00310 pccd_M->totface = cmd->numfaces; 00311 pccd_M->savety = CCD_SAVETY; 00312 pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f; 00313 pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f; 00314 pccd_M->mprevvert=NULL; 00315 00316 00317 /* blow it up with forcefield ranges */ 00318 hull = MAX2(ob->pd->pdef_sbift,ob->pd->pdef_sboft); 00319 00320 /* alloc and copy verts*/ 00321 pccd_M->mvert = MEM_dupallocN(cmd->xnew); 00322 /* note that xnew coords are already in global space, */ 00323 /* determine the ortho BB */ 00324 for(i=0; i < pccd_M->totvert; i++){ 00325 /* evaluate limits */ 00326 VECCOPY(v,pccd_M->mvert[i].co); 00327 pccd_M->bbmin[0] = MIN2(pccd_M->bbmin[0],v[0]-hull); 00328 pccd_M->bbmin[1] = MIN2(pccd_M->bbmin[1],v[1]-hull); 00329 pccd_M->bbmin[2] = MIN2(pccd_M->bbmin[2],v[2]-hull); 00330 00331 pccd_M->bbmax[0] = MAX2(pccd_M->bbmax[0],v[0]+hull); 00332 pccd_M->bbmax[1] = MAX2(pccd_M->bbmax[1],v[1]+hull); 00333 pccd_M->bbmax[2] = MAX2(pccd_M->bbmax[2],v[2]+hull); 00334 00335 } 00336 /* alloc and copy faces*/ 00337 pccd_M->mface = MEM_dupallocN(cmd->mfaces); 00338 00339 /* OBBs for idea1 */ 00340 pccd_M->mima = MEM_mallocN(sizeof(ccdf_minmax)*pccd_M->totface,"ccd_Mesh_Faces_mima"); 00341 mima = pccd_M->mima; 00342 mface = pccd_M->mface; 00343 00344 00345 /* anyhoo we need to walk the list of faces and find OBB they live in */ 00346 for(i=0; i < pccd_M->totface; i++){ 00347 mima->minx=mima->miny=mima->minz=1e30f; 00348 mima->maxx=mima->maxy=mima->maxz=-1e30f; 00349 00350 VECCOPY(v,pccd_M->mvert[mface->v1].co); 00351 mima->minx = MIN2(mima->minx,v[0]-hull); 00352 mima->miny = MIN2(mima->miny,v[1]-hull); 00353 mima->minz = MIN2(mima->minz,v[2]-hull); 00354 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00355 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00356 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00357 00358 VECCOPY(v,pccd_M->mvert[mface->v2].co); 00359 mima->minx = MIN2(mima->minx,v[0]-hull); 00360 mima->miny = MIN2(mima->miny,v[1]-hull); 00361 mima->minz = MIN2(mima->minz,v[2]-hull); 00362 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00363 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00364 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00365 00366 VECCOPY(v,pccd_M->mvert[mface->v3].co); 00367 mima->minx = MIN2(mima->minx,v[0]-hull); 00368 mima->miny = MIN2(mima->miny,v[1]-hull); 00369 mima->minz = MIN2(mima->minz,v[2]-hull); 00370 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00371 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00372 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00373 00374 if(mface->v4){ 00375 VECCOPY(v,pccd_M->mvert[mface->v4].co); 00376 mima->minx = MIN2(mima->minx,v[0]-hull); 00377 mima->miny = MIN2(mima->miny,v[1]-hull); 00378 mima->minz = MIN2(mima->minz,v[2]-hull); 00379 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00380 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00381 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00382 } 00383 00384 00385 mima++; 00386 mface++; 00387 00388 } 00389 return pccd_M; 00390 } 00391 static void ccd_mesh_update(Object *ob,ccd_Mesh *pccd_M) 00392 { 00393 CollisionModifierData *cmd; 00394 ccdf_minmax *mima =NULL; 00395 MFace *mface=NULL; 00396 float v[3],hull; 00397 int i; 00398 00399 cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision); 00400 00401 /* first some paranoia checks */ 00402 if (!cmd) return ; 00403 if (!cmd->numverts || !cmd->numfaces) return ; 00404 00405 if ((pccd_M->totvert != cmd->numverts) || 00406 (pccd_M->totface != cmd->numfaces)) return; 00407 00408 pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f; 00409 pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f; 00410 00411 00412 /* blow it up with forcefield ranges */ 00413 hull = MAX2(ob->pd->pdef_sbift,ob->pd->pdef_sboft); 00414 00415 /* rotate current to previous */ 00416 if(pccd_M->mprevvert) MEM_freeN(pccd_M->mprevvert); 00417 pccd_M->mprevvert = pccd_M->mvert; 00418 /* alloc and copy verts*/ 00419 pccd_M->mvert = MEM_dupallocN(cmd->xnew); 00420 /* note that xnew coords are already in global space, */ 00421 /* determine the ortho BB */ 00422 for(i=0; i < pccd_M->totvert; i++){ 00423 /* evaluate limits */ 00424 VECCOPY(v,pccd_M->mvert[i].co); 00425 pccd_M->bbmin[0] = MIN2(pccd_M->bbmin[0],v[0]-hull); 00426 pccd_M->bbmin[1] = MIN2(pccd_M->bbmin[1],v[1]-hull); 00427 pccd_M->bbmin[2] = MIN2(pccd_M->bbmin[2],v[2]-hull); 00428 00429 pccd_M->bbmax[0] = MAX2(pccd_M->bbmax[0],v[0]+hull); 00430 pccd_M->bbmax[1] = MAX2(pccd_M->bbmax[1],v[1]+hull); 00431 pccd_M->bbmax[2] = MAX2(pccd_M->bbmax[2],v[2]+hull); 00432 00433 /* evaluate limits */ 00434 VECCOPY(v,pccd_M->mprevvert[i].co); 00435 pccd_M->bbmin[0] = MIN2(pccd_M->bbmin[0],v[0]-hull); 00436 pccd_M->bbmin[1] = MIN2(pccd_M->bbmin[1],v[1]-hull); 00437 pccd_M->bbmin[2] = MIN2(pccd_M->bbmin[2],v[2]-hull); 00438 00439 pccd_M->bbmax[0] = MAX2(pccd_M->bbmax[0],v[0]+hull); 00440 pccd_M->bbmax[1] = MAX2(pccd_M->bbmax[1],v[1]+hull); 00441 pccd_M->bbmax[2] = MAX2(pccd_M->bbmax[2],v[2]+hull); 00442 00443 } 00444 00445 mima = pccd_M->mima; 00446 mface = pccd_M->mface; 00447 00448 00449 /* anyhoo we need to walk the list of faces and find OBB they live in */ 00450 for(i=0; i < pccd_M->totface; i++){ 00451 mima->minx=mima->miny=mima->minz=1e30f; 00452 mima->maxx=mima->maxy=mima->maxz=-1e30f; 00453 00454 VECCOPY(v,pccd_M->mvert[mface->v1].co); 00455 mima->minx = MIN2(mima->minx,v[0]-hull); 00456 mima->miny = MIN2(mima->miny,v[1]-hull); 00457 mima->minz = MIN2(mima->minz,v[2]-hull); 00458 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00459 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00460 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00461 00462 VECCOPY(v,pccd_M->mvert[mface->v2].co); 00463 mima->minx = MIN2(mima->minx,v[0]-hull); 00464 mima->miny = MIN2(mima->miny,v[1]-hull); 00465 mima->minz = MIN2(mima->minz,v[2]-hull); 00466 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00467 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00468 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00469 00470 VECCOPY(v,pccd_M->mvert[mface->v3].co); 00471 mima->minx = MIN2(mima->minx,v[0]-hull); 00472 mima->miny = MIN2(mima->miny,v[1]-hull); 00473 mima->minz = MIN2(mima->minz,v[2]-hull); 00474 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00475 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00476 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00477 00478 if(mface->v4){ 00479 VECCOPY(v,pccd_M->mvert[mface->v4].co); 00480 mima->minx = MIN2(mima->minx,v[0]-hull); 00481 mima->miny = MIN2(mima->miny,v[1]-hull); 00482 mima->minz = MIN2(mima->minz,v[2]-hull); 00483 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00484 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00485 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00486 } 00487 00488 00489 VECCOPY(v,pccd_M->mprevvert[mface->v1].co); 00490 mima->minx = MIN2(mima->minx,v[0]-hull); 00491 mima->miny = MIN2(mima->miny,v[1]-hull); 00492 mima->minz = MIN2(mima->minz,v[2]-hull); 00493 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00494 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00495 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00496 00497 VECCOPY(v,pccd_M->mprevvert[mface->v2].co); 00498 mima->minx = MIN2(mima->minx,v[0]-hull); 00499 mima->miny = MIN2(mima->miny,v[1]-hull); 00500 mima->minz = MIN2(mima->minz,v[2]-hull); 00501 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00502 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00503 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00504 00505 VECCOPY(v,pccd_M->mprevvert[mface->v3].co); 00506 mima->minx = MIN2(mima->minx,v[0]-hull); 00507 mima->miny = MIN2(mima->miny,v[1]-hull); 00508 mima->minz = MIN2(mima->minz,v[2]-hull); 00509 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00510 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00511 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00512 00513 if(mface->v4){ 00514 VECCOPY(v,pccd_M->mprevvert[mface->v4].co); 00515 mima->minx = MIN2(mima->minx,v[0]-hull); 00516 mima->miny = MIN2(mima->miny,v[1]-hull); 00517 mima->minz = MIN2(mima->minz,v[2]-hull); 00518 mima->maxx = MAX2(mima->maxx,v[0]+hull); 00519 mima->maxy = MAX2(mima->maxy,v[1]+hull); 00520 mima->maxz = MAX2(mima->maxz,v[2]+hull); 00521 } 00522 00523 00524 mima++; 00525 mface++; 00526 00527 } 00528 return ; 00529 } 00530 00531 static void ccd_mesh_free(ccd_Mesh *ccdm) 00532 { 00533 if(ccdm && (ccdm->savety == CCD_SAVETY )){ /*make sure we're not nuking objects we don't know*/ 00534 MEM_freeN(ccdm->mface); 00535 MEM_freeN(ccdm->mvert); 00536 if (ccdm->mprevvert) MEM_freeN(ccdm->mprevvert); 00537 MEM_freeN(ccdm->mima); 00538 MEM_freeN(ccdm); 00539 ccdm = NULL; 00540 } 00541 } 00542 00543 static void ccd_build_deflector_hash(Scene *scene, Object *vertexowner, GHash *hash) 00544 { 00545 Base *base= scene->base.first; 00546 Object *ob; 00547 00548 if (!hash) return; 00549 while (base) { 00550 /*Only proceed for mesh object in same layer */ 00551 if(base->object->type==OB_MESH && (base->lay & vertexowner->lay)) { 00552 ob= base->object; 00553 if((vertexowner) && (ob == vertexowner)) { 00554 /* if vertexowner is given we don't want to check collision with owner object */ 00555 base = base->next; 00556 continue; 00557 } 00558 00559 /*+++ only with deflecting set */ 00560 if(ob->pd && ob->pd->deflect && BLI_ghash_lookup(hash, ob) == NULL) { 00561 ccd_Mesh *ccdmesh = ccd_mesh_make(ob); 00562 BLI_ghash_insert(hash, ob, ccdmesh); 00563 }/*--- only with deflecting set */ 00564 00565 }/* mesh && layer*/ 00566 base = base->next; 00567 } /* while (base) */ 00568 } 00569 00570 static void ccd_update_deflector_hash(Scene *scene, Object *vertexowner, GHash *hash) 00571 { 00572 Base *base= scene->base.first; 00573 Object *ob; 00574 00575 if ((!hash) || (!vertexowner)) return; 00576 while (base) { 00577 /*Only proceed for mesh object in same layer */ 00578 if(base->object->type==OB_MESH && (base->lay & vertexowner->lay)) { 00579 ob= base->object; 00580 if(ob == vertexowner){ 00581 /* if vertexowner is given we don't want to check collision with owner object */ 00582 base = base->next; 00583 continue; 00584 } 00585 00586 /*+++ only with deflecting set */ 00587 if(ob->pd && ob->pd->deflect) { 00588 ccd_Mesh *ccdmesh = BLI_ghash_lookup(hash,ob); 00589 if (ccdmesh) 00590 ccd_mesh_update(ob,ccdmesh); 00591 }/*--- only with deflecting set */ 00592 00593 }/* mesh && layer*/ 00594 base = base->next; 00595 } /* while (base) */ 00596 } 00597 00598 00599 00600 00601 /*--- collider caching and dicing ---*/ 00602 00603 00604 static int count_mesh_quads(Mesh *me) 00605 { 00606 int a,result = 0; 00607 MFace *mface= me->mface; 00608 00609 if(mface) { 00610 for(a=me->totface; a>0; a--, mface++) { 00611 if(mface->v4) result++; 00612 } 00613 } 00614 return result; 00615 } 00616 00617 static void add_mesh_quad_diag_springs(Object *ob) 00618 { 00619 Mesh *me= ob->data; 00620 MFace *mface= me->mface; 00621 /*BodyPoint *bp;*/ /*UNUSED*/ 00622 BodySpring *bs, *bs_new; 00623 int a ; 00624 00625 if (ob->soft){ 00626 int nofquads; 00627 //float s_shear = ob->soft->shearstiff*ob->soft->shearstiff; 00628 00629 nofquads = count_mesh_quads(me); 00630 if (nofquads) { 00631 /* resize spring-array to hold additional quad springs */ 00632 bs_new= MEM_callocN( (ob->soft->totspring + nofquads *2 )*sizeof(BodySpring), "bodyspring"); 00633 memcpy(bs_new,ob->soft->bspring,(ob->soft->totspring )*sizeof(BodySpring)); 00634 00635 if(ob->soft->bspring) 00636 MEM_freeN(ob->soft->bspring); /* do this before reassigning the pointer or have a 1st class memory leak */ 00637 ob->soft->bspring = bs_new; 00638 00639 /* fill the tail */ 00640 a = 0; 00641 bs = bs_new+ob->soft->totspring; 00642 /*bp= ob->soft->bpoint; */ /*UNUSED*/ 00643 if(mface ) { 00644 for(a=me->totface; a>0; a--, mface++) { 00645 if(mface->v4) { 00646 bs->v1= mface->v1; 00647 bs->v2= mface->v3; 00648 bs->springtype =SB_STIFFQUAD; 00649 bs++; 00650 bs->v1= mface->v2; 00651 bs->v2= mface->v4; 00652 bs->springtype =SB_STIFFQUAD; 00653 bs++; 00654 00655 } 00656 } 00657 } 00658 00659 /* now we can announce new springs */ 00660 ob->soft->totspring += nofquads *2; 00661 } 00662 } 00663 } 00664 00665 static void add_2nd_order_roller(Object *ob,float UNUSED(stiffness), int *counter, int addsprings) 00666 { 00667 /*assume we have a softbody*/ 00668 SoftBody *sb= ob->soft; /* is supposed to be there */ 00669 BodyPoint *bp,*bpo; 00670 BodySpring *bs,*bs2,*bs3= NULL; 00671 int a,b,c,notthis= 0,v0; 00672 if (!sb->bspring){return;} /* we are 2nd order here so 1rst should have been build :) */ 00673 /* first run counting second run adding */ 00674 *counter = 0; 00675 if (addsprings) bs3 = ob->soft->bspring+ob->soft->totspring; 00676 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 00677 /*scan for neighborhood*/ 00678 bpo = NULL; 00679 v0 = (sb->totpoint-a); 00680 for(b=bp->nofsprings;b>0;b--){ 00681 bs = sb->bspring + bp->springs[b-1]; 00682 /*nasty thing here that springs have two ends 00683 so here we have to make sure we examine the other */ 00684 if (v0 == bs->v1){ 00685 bpo =sb->bpoint+bs->v2; 00686 notthis = bs->v2; 00687 } 00688 else { 00689 if (v0 == bs->v2){ 00690 bpo =sb->bpoint+bs->v1; 00691 notthis = bs->v1; 00692 } 00693 else {printf("oops we should not get here - add_2nd_order_springs");} 00694 } 00695 if (bpo){/* so now we have a 2nd order humpdidump */ 00696 for(c=bpo->nofsprings;c>0;c--){ 00697 bs2 = sb->bspring + bpo->springs[c-1]; 00698 if ((bs2->v1 != notthis) && (bs2->v1 > v0)){ 00699 (*counter)++;/*hit */ 00700 if (addsprings){ 00701 bs3->v1= v0; 00702 bs3->v2= bs2->v1; 00703 bs3->springtype =SB_BEND; 00704 bs3++; 00705 } 00706 } 00707 if ((bs2->v2 !=notthis)&&(bs2->v2 > v0)){ 00708 (*counter)++;/*hit */ 00709 if (addsprings){ 00710 bs3->v1= v0; 00711 bs3->v2= bs2->v2; 00712 bs3->springtype =SB_BEND; 00713 bs3++; 00714 } 00715 00716 } 00717 } 00718 00719 } 00720 00721 } 00722 /*scan for neighborhood done*/ 00723 } 00724 } 00725 00726 00727 static void add_2nd_order_springs(Object *ob,float stiffness) 00728 { 00729 int counter = 0; 00730 BodySpring *bs_new; 00731 stiffness *=stiffness; 00732 00733 add_2nd_order_roller(ob,stiffness,&counter,0); /* counting */ 00734 if (counter) { 00735 /* resize spring-array to hold additional springs */ 00736 bs_new= MEM_callocN( (ob->soft->totspring + counter )*sizeof(BodySpring), "bodyspring"); 00737 memcpy(bs_new,ob->soft->bspring,(ob->soft->totspring )*sizeof(BodySpring)); 00738 00739 if(ob->soft->bspring) 00740 MEM_freeN(ob->soft->bspring); 00741 ob->soft->bspring = bs_new; 00742 00743 add_2nd_order_roller(ob,stiffness,&counter,1); /* adding */ 00744 ob->soft->totspring +=counter ; 00745 } 00746 } 00747 00748 static void add_bp_springlist(BodyPoint *bp,int springID) 00749 { 00750 int *newlist; 00751 00752 if (bp->springs == NULL) { 00753 bp->springs = MEM_callocN( sizeof(int), "bpsprings"); 00754 bp->springs[0] = springID; 00755 bp->nofsprings = 1; 00756 } 00757 else { 00758 bp->nofsprings++; 00759 newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings"); 00760 memcpy(newlist,bp->springs,(bp->nofsprings-1)* sizeof(int)); 00761 MEM_freeN(bp->springs); 00762 bp->springs = newlist; 00763 bp->springs[bp->nofsprings-1] = springID; 00764 } 00765 } 00766 00767 /* do this once when sb is build 00768 it is O(N^2) so scanning for springs every iteration is too expensive 00769 */ 00770 static void build_bps_springlist(Object *ob) 00771 { 00772 SoftBody *sb= ob->soft; /* is supposed to be there */ 00773 BodyPoint *bp; 00774 BodySpring *bs; 00775 int a,b; 00776 00777 if (sb==NULL) return; /* paranoya check */ 00778 00779 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 00780 /* throw away old list */ 00781 if (bp->springs) { 00782 MEM_freeN(bp->springs); 00783 bp->springs=NULL; 00784 } 00785 /* scan for attached inner springs */ 00786 for(b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) { 00787 if (( (sb->totpoint-a) == bs->v1) ){ 00788 add_bp_springlist(bp,sb->totspring -b); 00789 } 00790 if (( (sb->totpoint-a) == bs->v2) ){ 00791 add_bp_springlist(bp,sb->totspring -b); 00792 } 00793 }/*for springs*/ 00794 }/*for bp*/ 00795 } 00796 00797 static void calculate_collision_balls(Object *ob) 00798 { 00799 SoftBody *sb= ob->soft; /* is supposed to be there */ 00800 BodyPoint *bp; 00801 BodySpring *bs; 00802 int a,b,akku_count; 00803 float min,max,akku; 00804 00805 if (sb==NULL) return; /* paranoya check */ 00806 00807 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 00808 bp->colball=0; 00809 akku =0.0f; 00810 akku_count=0; 00811 min = 1e22f; 00812 max = -1e22f; 00813 /* first estimation based on attached */ 00814 for(b=bp->nofsprings;b>0;b--){ 00815 bs = sb->bspring + bp->springs[b-1]; 00816 if (bs->springtype == SB_EDGE){ 00817 akku += bs->len; 00818 akku_count++, 00819 min = MIN2(bs->len,min); 00820 max = MAX2(bs->len,max); 00821 } 00822 } 00823 00824 if (akku_count > 0) { 00825 if (sb->sbc_mode == SBC_MODE_MANUAL){ 00826 bp->colball=sb->colball; 00827 } 00828 if (sb->sbc_mode == SBC_MODE_AVG){ 00829 bp->colball = akku/(float)akku_count*sb->colball; 00830 } 00831 if (sb->sbc_mode == SBC_MODE_MIN){ 00832 bp->colball=min*sb->colball; 00833 } 00834 if (sb->sbc_mode == SBC_MODE_MAX){ 00835 bp->colball=max*sb->colball; 00836 } 00837 if (sb->sbc_mode == SBC_MODE_AVGMINMAX){ 00838 bp->colball = (min + max)/2.0f*sb->colball; 00839 } 00840 } 00841 else bp->colball=0; 00842 }/*for bp*/ 00843 } 00844 00845 00846 /* creates new softbody if didn't exist yet, makes new points and springs arrays */ 00847 static void renew_softbody(Scene *scene, Object *ob, int totpoint, int totspring) 00848 { 00849 SoftBody *sb; 00850 int i; 00851 short softflag; 00852 if(ob->soft==NULL) ob->soft= sbNew(scene); 00853 else free_softbody_intern(ob->soft); 00854 sb= ob->soft; 00855 softflag=ob->softflag; 00856 00857 if(totpoint) { 00858 sb->totpoint= totpoint; 00859 sb->totspring= totspring; 00860 00861 sb->bpoint= MEM_mallocN( totpoint*sizeof(BodyPoint), "bodypoint"); 00862 if(totspring) 00863 sb->bspring= MEM_mallocN( totspring*sizeof(BodySpring), "bodyspring"); 00864 00865 /* initialise BodyPoint array */ 00866 for (i=0; i<totpoint; i++) { 00867 BodyPoint *bp = &sb->bpoint[i]; 00868 00869 00870 /* hum as far as i see this is overridden by _final_goal() now jow_go_for2_5 */ 00871 /* sadly breaks compatibility with older versions */ 00872 /* but makes goals behave the same for meshes, lattices and curves */ 00873 if(softflag & OB_SB_GOAL) { 00874 bp->goal= sb->defgoal; 00875 } 00876 else { 00877 bp->goal= 0.0f; 00878 /* so this will definily be below SOFTGOALSNAP */ 00879 } 00880 00881 bp->nofsprings= 0; 00882 bp->springs= NULL; 00883 bp->choke = 0.0f; 00884 bp->choke2 = 0.0f; 00885 bp->frozen = 1.0f; 00886 bp->colball = 0.0f; 00887 bp->loc_flag = 0; 00888 bp->springweight = 1.0f; 00889 bp->mass = 1.0f; 00890 } 00891 } 00892 } 00893 00894 static void free_softbody_baked(SoftBody *sb) 00895 { 00896 SBVertex *key; 00897 int k; 00898 00899 for(k=0; k<sb->totkey; k++) { 00900 key= *(sb->keys + k); 00901 if(key) MEM_freeN(key); 00902 } 00903 if(sb->keys) MEM_freeN(sb->keys); 00904 00905 sb->keys= NULL; 00906 sb->totkey= 0; 00907 } 00908 static void free_scratch(SoftBody *sb) 00909 { 00910 if(sb->scratch){ 00911 /* todo make sure everything is cleaned up nicly */ 00912 if (sb->scratch->colliderhash){ 00913 BLI_ghash_free(sb->scratch->colliderhash, NULL, 00914 (GHashValFreeFP) ccd_mesh_free); /*this hoepfully will free all caches*/ 00915 sb->scratch->colliderhash = NULL; 00916 } 00917 if (sb->scratch->bodyface){ 00918 MEM_freeN(sb->scratch->bodyface); 00919 } 00920 if (sb->scratch->Ref.ivert){ 00921 MEM_freeN(sb->scratch->Ref.ivert); 00922 } 00923 MEM_freeN(sb->scratch); 00924 sb->scratch = NULL; 00925 } 00926 00927 } 00928 00929 /* only frees internal data */ 00930 static void free_softbody_intern(SoftBody *sb) 00931 { 00932 if(sb) { 00933 int a; 00934 BodyPoint *bp; 00935 00936 if(sb->bpoint){ 00937 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 00938 /* free spring list */ 00939 if (bp->springs != NULL) { 00940 MEM_freeN(bp->springs); 00941 } 00942 } 00943 MEM_freeN(sb->bpoint); 00944 } 00945 00946 if(sb->bspring) MEM_freeN(sb->bspring); 00947 00948 sb->totpoint= sb->totspring= 0; 00949 sb->bpoint= NULL; 00950 sb->bspring= NULL; 00951 00952 free_scratch(sb); 00953 free_softbody_baked(sb); 00954 } 00955 } 00956 00957 00958 /* ************ dynamics ********** */ 00959 00960 /* the most general (micro physics correct) way to do collision 00961 ** (only needs the current particle position) 00962 ** 00963 ** it actually checks if the particle intrudes a short range force field generated 00964 ** by the faces of the target object and returns a force to drive the particel out 00965 ** the strenght of the field grows exponetially if the particle is on the 'wrong' side of the face 00966 ** 'wrong' side : projection to the face normal is negative (all referred to a vertex in the face) 00967 ** 00968 ** flaw of this: 'fast' particles as well as 'fast' colliding faces 00969 ** give a 'tunnel' effect such that the particle passes through the force field 00970 ** without ever 'seeing' it 00971 ** this is fully compliant to heisenberg: h >= fuzzy(location) * fuzzy(time) 00972 ** besides our h is way larger than in QM because forces propagate way slower here 00973 ** we have to deal with fuzzy(time) in the range of 1/25 seconds (typical frame rate) 00974 ** yup collision targets are not known here any better 00975 ** and 1/25 second is looong compared to real collision events 00976 ** Q: why not use 'simple' collision here like bouncing back a particle 00977 ** --> reverting is velocity on the face normal 00978 ** A: because our particles are not alone here 00979 ** and need to tell their neighbours exactly what happens via spring forces 00980 ** unless sbObjectStep( .. ) is called on sub frame timing level 00981 ** BTW that also questions the use of a 'implicit' solvers on softbodies 00982 ** since that would only valid for 'slow' moving collision targets and dito particles 00983 */ 00984 00985 /* aye this belongs to arith.c */ 00986 static void Vec3PlusStVec(float *v, float s, float *v1) 00987 { 00988 v[0] += s*v1[0]; 00989 v[1] += s*v1[1]; 00990 v[2] += s*v1[2]; 00991 } 00992 00993 /* +++ dependancy information functions*/ 00994 00995 static int are_there_deflectors(Scene *scene, unsigned int layer) 00996 { 00997 Base *base; 00998 00999 for(base = scene->base.first; base; base= base->next) { 01000 if( (base->lay & layer) && base->object->pd) { 01001 if(base->object->pd->deflect) 01002 return 1; 01003 } 01004 } 01005 return 0; 01006 } 01007 01008 static int query_external_colliders(Scene *scene, Object *me) 01009 { 01010 return(are_there_deflectors(scene, me->lay)); 01011 } 01012 /* --- dependancy information functions*/ 01013 01014 01015 /* +++ the aabb "force" section*/ 01016 static int sb_detect_aabb_collisionCached( float UNUSED(force[3]), unsigned int UNUSED(par_layer),struct Object *vertexowner,float UNUSED(time)) 01017 { 01018 Object *ob; 01019 SoftBody *sb=vertexowner->soft; 01020 GHash *hash; 01021 GHashIterator *ihash; 01022 float aabbmin[3],aabbmax[3]; 01023 int deflected=0; 01024 #if 0 01025 int a; 01026 #endif 01027 01028 if ((sb == NULL) || (sb->scratch ==NULL)) return 0; 01029 VECCOPY(aabbmin,sb->scratch->aabbmin); 01030 VECCOPY(aabbmax,sb->scratch->aabbmax); 01031 01032 hash = vertexowner->soft->scratch->colliderhash; 01033 ihash = BLI_ghashIterator_new(hash); 01034 while (!BLI_ghashIterator_isDone(ihash) ) { 01035 01036 ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); 01037 ob = BLI_ghashIterator_getKey (ihash); 01038 /* only with deflecting set */ 01039 if(ob->pd && ob->pd->deflect) { 01040 #if 0 /* UNUSED */ 01041 MFace *mface= NULL; 01042 MVert *mvert= NULL; 01043 MVert *mprevvert= NULL; 01044 ccdf_minmax *mima= NULL; 01045 #endif 01046 if(ccdm){ 01047 #if 0 /* UNUSED */ 01048 mface= ccdm->mface; 01049 mvert= ccdm->mvert; 01050 mprevvert= ccdm->mprevvert; 01051 mima= ccdm->mima; 01052 a = ccdm->totface; 01053 #endif 01054 if ((aabbmax[0] < ccdm->bbmin[0]) || 01055 (aabbmax[1] < ccdm->bbmin[1]) || 01056 (aabbmax[2] < ccdm->bbmin[2]) || 01057 (aabbmin[0] > ccdm->bbmax[0]) || 01058 (aabbmin[1] > ccdm->bbmax[1]) || 01059 (aabbmin[2] > ccdm->bbmax[2]) ) { 01060 /* boxes dont intersect */ 01061 BLI_ghashIterator_step(ihash); 01062 continue; 01063 } 01064 01065 /* so now we have the 2 boxes overlapping */ 01066 /* forces actually not used */ 01067 deflected = 2; 01068 01069 } 01070 else{ 01071 /*aye that should be cached*/ 01072 printf("missing cache error \n"); 01073 BLI_ghashIterator_step(ihash); 01074 continue; 01075 } 01076 } /* if(ob->pd && ob->pd->deflect) */ 01077 BLI_ghashIterator_step(ihash); 01078 } /* while () */ 01079 BLI_ghashIterator_free(ihash); 01080 return deflected; 01081 } 01082 /* --- the aabb section*/ 01083 01084 01085 /* +++ the face external section*/ 01086 static int sb_detect_face_pointCached(float face_v1[3],float face_v2[3],float face_v3[3],float *damp, 01087 float force[3], unsigned int UNUSED(par_layer),struct Object *vertexowner,float time) 01088 { 01089 Object *ob; 01090 GHash *hash; 01091 GHashIterator *ihash; 01092 float nv1[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3],aabbmax[3]; 01093 float facedist,outerfacethickness,tune = 10.f; 01094 int a, deflected=0; 01095 01096 aabbmin[0] = MIN3(face_v1[0],face_v2[0],face_v3[0]); 01097 aabbmin[1] = MIN3(face_v1[1],face_v2[1],face_v3[1]); 01098 aabbmin[2] = MIN3(face_v1[2],face_v2[2],face_v3[2]); 01099 aabbmax[0] = MAX3(face_v1[0],face_v2[0],face_v3[0]); 01100 aabbmax[1] = MAX3(face_v1[1],face_v2[1],face_v3[1]); 01101 aabbmax[2] = MAX3(face_v1[2],face_v2[2],face_v3[2]); 01102 01103 /* calculate face normal once again SIGH */ 01104 VECSUB(edge1, face_v1, face_v2); 01105 VECSUB(edge2, face_v3, face_v2); 01106 cross_v3_v3v3(d_nvect, edge2, edge1); 01107 normalize_v3(d_nvect); 01108 01109 01110 hash = vertexowner->soft->scratch->colliderhash; 01111 ihash = BLI_ghashIterator_new(hash); 01112 while (!BLI_ghashIterator_isDone(ihash) ) { 01113 01114 ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); 01115 ob = BLI_ghashIterator_getKey (ihash); 01116 /* only with deflecting set */ 01117 if(ob->pd && ob->pd->deflect) { 01118 MVert *mvert= NULL; 01119 MVert *mprevvert= NULL; 01120 if(ccdm){ 01121 mvert= ccdm->mvert; 01122 a = ccdm->totvert; 01123 mprevvert= ccdm->mprevvert; 01124 outerfacethickness =ob->pd->pdef_sboft; 01125 if ((aabbmax[0] < ccdm->bbmin[0]) || 01126 (aabbmax[1] < ccdm->bbmin[1]) || 01127 (aabbmax[2] < ccdm->bbmin[2]) || 01128 (aabbmin[0] > ccdm->bbmax[0]) || 01129 (aabbmin[1] > ccdm->bbmax[1]) || 01130 (aabbmin[2] > ccdm->bbmax[2]) ) { 01131 /* boxes dont intersect */ 01132 BLI_ghashIterator_step(ihash); 01133 continue; 01134 } 01135 01136 } 01137 else{ 01138 /*aye that should be cached*/ 01139 printf("missing cache error \n"); 01140 BLI_ghashIterator_step(ihash); 01141 continue; 01142 } 01143 01144 01145 /* use mesh*/ 01146 if (mvert) { 01147 while(a){ 01148 VECCOPY(nv1,mvert[a-1].co); 01149 if(mprevvert){ 01150 mul_v3_fl(nv1,time); 01151 Vec3PlusStVec(nv1,(1.0f-time),mprevvert[a-1].co); 01152 } 01153 /* origin to face_v2*/ 01154 VECSUB(nv1, nv1, face_v2); 01155 facedist = dot_v3v3(nv1,d_nvect); 01156 if (ABS(facedist)<outerfacethickness){ 01157 if (isect_point_tri_prism_v3(nv1, face_v1,face_v2,face_v3) ){ 01158 float df; 01159 if (facedist > 0){ 01160 df = (outerfacethickness-facedist)/outerfacethickness; 01161 } 01162 else { 01163 df = (outerfacethickness+facedist)/outerfacethickness; 01164 } 01165 01166 *damp=df*tune*ob->pd->pdef_sbdamp; 01167 01168 df = 0.01f*exp(- 100.0f*df); 01169 Vec3PlusStVec(force,-df,d_nvect); 01170 deflected = 3; 01171 } 01172 } 01173 a--; 01174 }/* while(a)*/ 01175 } /* if (mvert) */ 01176 } /* if(ob->pd && ob->pd->deflect) */ 01177 BLI_ghashIterator_step(ihash); 01178 } /* while () */ 01179 BLI_ghashIterator_free(ihash); 01180 return deflected; 01181 } 01182 01183 01184 static int sb_detect_face_collisionCached(float face_v1[3],float face_v2[3],float face_v3[3],float *damp, 01185 float force[3], unsigned int UNUSED(par_layer),struct Object *vertexowner,float time) 01186 { 01187 Object *ob; 01188 GHash *hash; 01189 GHashIterator *ihash; 01190 float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3],aabbmax[3]; 01191 float t,tune = 10.0f; 01192 int a, deflected=0; 01193 01194 aabbmin[0] = MIN3(face_v1[0],face_v2[0],face_v3[0]); 01195 aabbmin[1] = MIN3(face_v1[1],face_v2[1],face_v3[1]); 01196 aabbmin[2] = MIN3(face_v1[2],face_v2[2],face_v3[2]); 01197 aabbmax[0] = MAX3(face_v1[0],face_v2[0],face_v3[0]); 01198 aabbmax[1] = MAX3(face_v1[1],face_v2[1],face_v3[1]); 01199 aabbmax[2] = MAX3(face_v1[2],face_v2[2],face_v3[2]); 01200 01201 hash = vertexowner->soft->scratch->colliderhash; 01202 ihash = BLI_ghashIterator_new(hash); 01203 while (!BLI_ghashIterator_isDone(ihash) ) { 01204 01205 ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); 01206 ob = BLI_ghashIterator_getKey (ihash); 01207 /* only with deflecting set */ 01208 if(ob->pd && ob->pd->deflect) { 01209 MFace *mface= NULL; 01210 MVert *mvert= NULL; 01211 MVert *mprevvert= NULL; 01212 ccdf_minmax *mima= NULL; 01213 if(ccdm){ 01214 mface= ccdm->mface; 01215 mvert= ccdm->mvert; 01216 mprevvert= ccdm->mprevvert; 01217 mima= ccdm->mima; 01218 a = ccdm->totface; 01219 01220 if ((aabbmax[0] < ccdm->bbmin[0]) || 01221 (aabbmax[1] < ccdm->bbmin[1]) || 01222 (aabbmax[2] < ccdm->bbmin[2]) || 01223 (aabbmin[0] > ccdm->bbmax[0]) || 01224 (aabbmin[1] > ccdm->bbmax[1]) || 01225 (aabbmin[2] > ccdm->bbmax[2]) ) { 01226 /* boxes dont intersect */ 01227 BLI_ghashIterator_step(ihash); 01228 continue; 01229 } 01230 01231 } 01232 else{ 01233 /*aye that should be cached*/ 01234 printf("missing cache error \n"); 01235 BLI_ghashIterator_step(ihash); 01236 continue; 01237 } 01238 01239 01240 /* use mesh*/ 01241 while (a--) { 01242 if ( 01243 (aabbmax[0] < mima->minx) || 01244 (aabbmin[0] > mima->maxx) || 01245 (aabbmax[1] < mima->miny) || 01246 (aabbmin[1] > mima->maxy) || 01247 (aabbmax[2] < mima->minz) || 01248 (aabbmin[2] > mima->maxz) 01249 ) { 01250 mface++; 01251 mima++; 01252 continue; 01253 } 01254 01255 01256 if (mvert){ 01257 01258 VECCOPY(nv1,mvert[mface->v1].co); 01259 VECCOPY(nv2,mvert[mface->v2].co); 01260 VECCOPY(nv3,mvert[mface->v3].co); 01261 if (mface->v4){ 01262 VECCOPY(nv4,mvert[mface->v4].co); 01263 } 01264 if (mprevvert){ 01265 mul_v3_fl(nv1,time); 01266 Vec3PlusStVec(nv1,(1.0f-time),mprevvert[mface->v1].co); 01267 01268 mul_v3_fl(nv2,time); 01269 Vec3PlusStVec(nv2,(1.0f-time),mprevvert[mface->v2].co); 01270 01271 mul_v3_fl(nv3,time); 01272 Vec3PlusStVec(nv3,(1.0f-time),mprevvert[mface->v3].co); 01273 01274 if (mface->v4){ 01275 mul_v3_fl(nv4,time); 01276 Vec3PlusStVec(nv4,(1.0f-time),mprevvert[mface->v4].co); 01277 } 01278 } 01279 } 01280 01281 /* switch origin to be nv2*/ 01282 VECSUB(edge1, nv1, nv2); 01283 VECSUB(edge2, nv3, nv2); 01284 cross_v3_v3v3(d_nvect, edge2, edge1); 01285 normalize_v3(d_nvect); 01286 if ( 01287 isect_line_tri_v3(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) || 01288 isect_line_tri_v3(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) || 01289 isect_line_tri_v3(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL) ){ 01290 Vec3PlusStVec(force,-0.5f,d_nvect); 01291 *damp=tune*ob->pd->pdef_sbdamp; 01292 deflected = 2; 01293 } 01294 if (mface->v4){ /* quad */ 01295 /* switch origin to be nv4 */ 01296 VECSUB(edge1, nv3, nv4); 01297 VECSUB(edge2, nv1, nv4); 01298 cross_v3_v3v3(d_nvect, edge2, edge1); 01299 normalize_v3(d_nvect); 01300 if ( 01301 /* isect_line_tri_v3(nv1, nv3, face_v1, face_v2, face_v3, &t, NULL) || 01302 we did that edge already */ 01303 isect_line_tri_v3(nv3, nv4, face_v1, face_v2, face_v3, &t, NULL) || 01304 isect_line_tri_v3(nv4, nv1, face_v1, face_v2, face_v3, &t, NULL) ){ 01305 Vec3PlusStVec(force,-0.5f,d_nvect); 01306 *damp=tune*ob->pd->pdef_sbdamp; 01307 deflected = 2; 01308 } 01309 } 01310 mface++; 01311 mima++; 01312 }/* while a */ 01313 } /* if(ob->pd && ob->pd->deflect) */ 01314 BLI_ghashIterator_step(ihash); 01315 } /* while () */ 01316 BLI_ghashIterator_free(ihash); 01317 return deflected; 01318 } 01319 01320 01321 01322 static void scan_for_ext_face_forces(Object *ob,float timenow) 01323 { 01324 SoftBody *sb = ob->soft; 01325 BodyFace *bf; 01326 int a; 01327 float damp=0.0f,choke=1.0f; 01328 float tune = -10.0f; 01329 float feedback[3]; 01330 01331 if (sb && sb->scratch->totface){ 01332 01333 01334 bf = sb->scratch->bodyface; 01335 for(a=0; a<sb->scratch->totface; a++, bf++) { 01336 bf->ext_force[0]=bf->ext_force[1]=bf->ext_force[2]=0.0f; 01337 /*+++edges intruding*/ 01338 bf->flag &= ~BFF_INTERSECT; 01339 feedback[0]=feedback[1]=feedback[2]=0.0f; 01340 if (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, 01341 &damp, feedback, ob->lay ,ob , timenow)){ 01342 Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); 01343 Vec3PlusStVec(sb->bpoint[bf->v2].force,tune,feedback); 01344 Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); 01345 // Vec3PlusStVec(bf->ext_force,tune,feedback); 01346 bf->flag |= BFF_INTERSECT; 01347 choke = MIN2(MAX2(damp,choke),1.0f); 01348 } 01349 01350 feedback[0]=feedback[1]=feedback[2]=0.0f; 01351 if ((bf->v4) && (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, 01352 &damp, feedback, ob->lay ,ob , timenow))){ 01353 Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); 01354 Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); 01355 Vec3PlusStVec(sb->bpoint[bf->v4].force,tune,feedback); 01356 // Vec3PlusStVec(bf->ext_force,tune,feedback); 01357 bf->flag |= BFF_INTERSECT; 01358 choke = MIN2(MAX2(damp,choke),1.0f); 01359 } 01360 /*---edges intruding*/ 01361 01362 /*+++ close vertices*/ 01363 if (( bf->flag & BFF_INTERSECT)==0){ 01364 bf->flag &= ~BFF_CLOSEVERT; 01365 tune = -1.0f; 01366 feedback[0]=feedback[1]=feedback[2]=0.0f; 01367 if (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, 01368 &damp, feedback, ob->lay ,ob , timenow)){ 01369 Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); 01370 Vec3PlusStVec(sb->bpoint[bf->v2].force,tune,feedback); 01371 Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); 01372 // Vec3PlusStVec(bf->ext_force,tune,feedback); 01373 bf->flag |= BFF_CLOSEVERT; 01374 choke = MIN2(MAX2(damp,choke),1.0f); 01375 } 01376 01377 feedback[0]=feedback[1]=feedback[2]=0.0f; 01378 if ((bf->v4) && (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos,sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, 01379 &damp, feedback, ob->lay ,ob , timenow))){ 01380 Vec3PlusStVec(sb->bpoint[bf->v1].force,tune,feedback); 01381 Vec3PlusStVec(sb->bpoint[bf->v3].force,tune,feedback); 01382 Vec3PlusStVec(sb->bpoint[bf->v4].force,tune,feedback); 01383 // Vec3PlusStVec(bf->ext_force,tune,feedback); 01384 bf->flag |= BFF_CLOSEVERT; 01385 choke = MIN2(MAX2(damp,choke),1.0f); 01386 } 01387 } 01388 /*--- close vertices*/ 01389 } 01390 bf = sb->scratch->bodyface; 01391 for(a=0; a<sb->scratch->totface; a++, bf++) { 01392 if (( bf->flag & BFF_INTERSECT) || ( bf->flag & BFF_CLOSEVERT)) 01393 { 01394 sb->bpoint[bf->v1].choke2=MAX2(sb->bpoint[bf->v1].choke2,choke); 01395 sb->bpoint[bf->v2].choke2=MAX2(sb->bpoint[bf->v2].choke2,choke); 01396 sb->bpoint[bf->v3].choke2=MAX2(sb->bpoint[bf->v3].choke2,choke); 01397 if (bf->v4){ 01398 sb->bpoint[bf->v2].choke2=MAX2(sb->bpoint[bf->v2].choke2,choke); 01399 } 01400 } 01401 } 01402 } 01403 } 01404 01405 /* --- the face external section*/ 01406 01407 01408 /* +++ the spring external section*/ 01409 01410 static int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp, 01411 float force[3], unsigned int UNUSED(par_layer),struct Object *vertexowner,float time) 01412 { 01413 Object *ob; 01414 GHash *hash; 01415 GHashIterator *ihash; 01416 float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3],aabbmax[3]; 01417 float t,el; 01418 int a, deflected=0; 01419 01420 aabbmin[0] = MIN2(edge_v1[0],edge_v2[0]); 01421 aabbmin[1] = MIN2(edge_v1[1],edge_v2[1]); 01422 aabbmin[2] = MIN2(edge_v1[2],edge_v2[2]); 01423 aabbmax[0] = MAX2(edge_v1[0],edge_v2[0]); 01424 aabbmax[1] = MAX2(edge_v1[1],edge_v2[1]); 01425 aabbmax[2] = MAX2(edge_v1[2],edge_v2[2]); 01426 01427 el = len_v3v3(edge_v1,edge_v2); 01428 01429 hash = vertexowner->soft->scratch->colliderhash; 01430 ihash = BLI_ghashIterator_new(hash); 01431 while (!BLI_ghashIterator_isDone(ihash) ) { 01432 01433 ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); 01434 ob = BLI_ghashIterator_getKey (ihash); 01435 /* only with deflecting set */ 01436 if(ob->pd && ob->pd->deflect) { 01437 MFace *mface= NULL; 01438 MVert *mvert= NULL; 01439 MVert *mprevvert= NULL; 01440 ccdf_minmax *mima= NULL; 01441 if(ccdm){ 01442 mface= ccdm->mface; 01443 mvert= ccdm->mvert; 01444 mprevvert= ccdm->mprevvert; 01445 mima= ccdm->mima; 01446 a = ccdm->totface; 01447 01448 if ((aabbmax[0] < ccdm->bbmin[0]) || 01449 (aabbmax[1] < ccdm->bbmin[1]) || 01450 (aabbmax[2] < ccdm->bbmin[2]) || 01451 (aabbmin[0] > ccdm->bbmax[0]) || 01452 (aabbmin[1] > ccdm->bbmax[1]) || 01453 (aabbmin[2] > ccdm->bbmax[2]) ) { 01454 /* boxes dont intersect */ 01455 BLI_ghashIterator_step(ihash); 01456 continue; 01457 } 01458 01459 } 01460 else{ 01461 /*aye that should be cached*/ 01462 printf("missing cache error \n"); 01463 BLI_ghashIterator_step(ihash); 01464 continue; 01465 } 01466 01467 01468 /* use mesh*/ 01469 while (a--) { 01470 if ( 01471 (aabbmax[0] < mima->minx) || 01472 (aabbmin[0] > mima->maxx) || 01473 (aabbmax[1] < mima->miny) || 01474 (aabbmin[1] > mima->maxy) || 01475 (aabbmax[2] < mima->minz) || 01476 (aabbmin[2] > mima->maxz) 01477 ) { 01478 mface++; 01479 mima++; 01480 continue; 01481 } 01482 01483 01484 if (mvert){ 01485 01486 VECCOPY(nv1,mvert[mface->v1].co); 01487 VECCOPY(nv2,mvert[mface->v2].co); 01488 VECCOPY(nv3,mvert[mface->v3].co); 01489 if (mface->v4){ 01490 VECCOPY(nv4,mvert[mface->v4].co); 01491 } 01492 if (mprevvert){ 01493 mul_v3_fl(nv1,time); 01494 Vec3PlusStVec(nv1,(1.0f-time),mprevvert[mface->v1].co); 01495 01496 mul_v3_fl(nv2,time); 01497 Vec3PlusStVec(nv2,(1.0f-time),mprevvert[mface->v2].co); 01498 01499 mul_v3_fl(nv3,time); 01500 Vec3PlusStVec(nv3,(1.0f-time),mprevvert[mface->v3].co); 01501 01502 if (mface->v4){ 01503 mul_v3_fl(nv4,time); 01504 Vec3PlusStVec(nv4,(1.0f-time),mprevvert[mface->v4].co); 01505 } 01506 } 01507 } 01508 01509 /* switch origin to be nv2*/ 01510 VECSUB(edge1, nv1, nv2); 01511 VECSUB(edge2, nv3, nv2); 01512 01513 cross_v3_v3v3(d_nvect, edge2, edge1); 01514 normalize_v3(d_nvect); 01515 if ( isect_line_tri_v3(edge_v1, edge_v2, nv1, nv2, nv3, &t, NULL)){ 01516 float v1[3],v2[3]; 01517 float intrusiondepth,i1,i2; 01518 VECSUB(v1, edge_v1, nv2); 01519 VECSUB(v2, edge_v2, nv2); 01520 i1 = dot_v3v3(v1,d_nvect); 01521 i2 = dot_v3v3(v2,d_nvect); 01522 intrusiondepth = -MIN2(i1,i2)/el; 01523 Vec3PlusStVec(force,intrusiondepth,d_nvect); 01524 *damp=ob->pd->pdef_sbdamp; 01525 deflected = 2; 01526 } 01527 if (mface->v4){ /* quad */ 01528 /* switch origin to be nv4 */ 01529 VECSUB(edge1, nv3, nv4); 01530 VECSUB(edge2, nv1, nv4); 01531 01532 cross_v3_v3v3(d_nvect, edge2, edge1); 01533 normalize_v3(d_nvect); 01534 if (isect_line_tri_v3( edge_v1, edge_v2,nv1, nv3, nv4, &t, NULL)){ 01535 float v1[3],v2[3]; 01536 float intrusiondepth,i1,i2; 01537 VECSUB(v1, edge_v1, nv4); 01538 VECSUB(v2, edge_v2, nv4); 01539 i1 = dot_v3v3(v1,d_nvect); 01540 i2 = dot_v3v3(v2,d_nvect); 01541 intrusiondepth = -MIN2(i1,i2)/el; 01542 01543 01544 Vec3PlusStVec(force,intrusiondepth,d_nvect); 01545 *damp=ob->pd->pdef_sbdamp; 01546 deflected = 2; 01547 } 01548 } 01549 mface++; 01550 mima++; 01551 }/* while a */ 01552 } /* if(ob->pd && ob->pd->deflect) */ 01553 BLI_ghashIterator_step(ihash); 01554 } /* while () */ 01555 BLI_ghashIterator_free(ihash); 01556 return deflected; 01557 } 01558 01559 static void _scan_for_ext_spring_forces(Scene *scene, Object *ob, float timenow, int ifirst, int ilast, struct ListBase *do_effector) 01560 { 01561 SoftBody *sb = ob->soft; 01562 int a; 01563 float damp; 01564 float feedback[3]; 01565 01566 if (sb && sb->totspring){ 01567 for(a=ifirst; a<ilast; a++) { 01568 BodySpring *bs = &sb->bspring[a]; 01569 bs->ext_force[0]=bs->ext_force[1]=bs->ext_force[2]=0.0f; 01570 feedback[0]=feedback[1]=feedback[2]=0.0f; 01571 bs->flag &= ~BSF_INTERSECT; 01572 01573 if (bs->springtype == SB_EDGE){ 01574 /* +++ springs colliding */ 01575 if (ob->softflag & OB_SB_EDGECOLL){ 01576 if ( sb_detect_edge_collisionCached (sb->bpoint[bs->v1].pos , sb->bpoint[bs->v2].pos, 01577 &damp,feedback,ob->lay,ob,timenow)){ 01578 add_v3_v3(bs->ext_force, feedback); 01579 bs->flag |= BSF_INTERSECT; 01580 //bs->cf=damp; 01581 bs->cf=sb->choke*0.01f; 01582 01583 } 01584 } 01585 /* ---- springs colliding */ 01586 01587 /* +++ springs seeing wind ... n stuff depending on their orientation*/ 01588 /* note we don't use sb->mediafrict but use sb->aeroedge for magnitude of effect*/ 01589 if(sb->aeroedge){ 01590 float vel[3],sp[3],pr[3],force[3]; 01591 float f,windfactor = 0.25f; 01592 /*see if we have wind*/ 01593 if(do_effector) { 01594 EffectedPoint epoint; 01595 float speed[3]={0.0f,0.0f,0.0f}; 01596 float pos[3]; 01597 mid_v3_v3v3(pos, sb->bpoint[bs->v1].pos , sb->bpoint[bs->v2].pos); 01598 mid_v3_v3v3(vel, sb->bpoint[bs->v1].vec , sb->bpoint[bs->v2].vec); 01599 pd_point_from_soft(scene, pos, vel, -1, &epoint); 01600 pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); 01601 01602 mul_v3_fl(speed,windfactor); 01603 add_v3_v3(vel, speed); 01604 } 01605 /* media in rest */ 01606 else{ 01607 VECADD(vel, sb->bpoint[bs->v1].vec , sb->bpoint[bs->v2].vec); 01608 } 01609 f = normalize_v3(vel); 01610 f = -0.0001f*f*f*sb->aeroedge; 01611 /* (todo) add a nice angle dependant function done for now BUT */ 01612 /* still there could be some nice drag/lift function, but who needs it */ 01613 01614 VECSUB(sp, sb->bpoint[bs->v1].pos , sb->bpoint[bs->v2].pos); 01615 project_v3_v3v3(pr,vel,sp); 01616 VECSUB(vel,vel,pr); 01617 normalize_v3(vel); 01618 if (ob->softflag & OB_SB_AERO_ANGLE){ 01619 normalize_v3(sp); 01620 Vec3PlusStVec(bs->ext_force,f*(1.0f-ABS(dot_v3v3(vel,sp))),vel); 01621 } 01622 else{ 01623 Vec3PlusStVec(bs->ext_force,f,vel); // to keep compatible with 2.45 release files 01624 } 01625 } 01626 /* --- springs seeing wind */ 01627 } 01628 } 01629 } 01630 } 01631 01632 01633 static void scan_for_ext_spring_forces(Scene *scene, Object *ob, float timenow) 01634 { 01635 SoftBody *sb = ob->soft; 01636 ListBase *do_effector = NULL; 01637 01638 do_effector = pdInitEffectors(scene, ob, NULL, sb->effector_weights); 01639 _scan_for_ext_spring_forces(scene, ob, timenow, 0, sb->totspring, do_effector); 01640 pdEndEffectors(&do_effector); 01641 } 01642 01643 static void *exec_scan_for_ext_spring_forces(void *data) 01644 { 01645 SB_thread_context *pctx = (SB_thread_context*)data; 01646 _scan_for_ext_spring_forces(pctx->scene, pctx->ob, pctx->timenow, pctx->ifirst, pctx->ilast, pctx->do_effector); 01647 return NULL; 01648 } 01649 01650 static void sb_sfesf_threads_run(Scene *scene, struct Object *ob, float timenow,int totsprings,int *UNUSED(ptr_to_break_func(void))) 01651 { 01652 ListBase *do_effector = NULL; 01653 ListBase threads; 01654 SB_thread_context *sb_threads; 01655 int i, totthread,left,dec; 01656 int lowsprings =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ 01657 01658 do_effector= pdInitEffectors(scene, ob, NULL, ob->soft->effector_weights); 01659 01660 /* figure the number of threads while preventing pretty pointless threading overhead */ 01661 if(scene->r.mode & R_FIXED_THREADS) 01662 totthread= scene->r.threads; 01663 else 01664 totthread= BLI_system_thread_count(); 01665 /* what if we got zillions of CPUs running but less to spread*/ 01666 while ((totsprings/totthread < lowsprings) && (totthread > 1)) { 01667 totthread--; 01668 } 01669 01670 sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBSpringsThread"); 01671 memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); 01672 left = totsprings; 01673 dec = totsprings/totthread +1; 01674 for(i=0; i<totthread; i++) { 01675 sb_threads[i].scene = scene; 01676 sb_threads[i].ob = ob; 01677 sb_threads[i].forcetime = 0.0; // not used here 01678 sb_threads[i].timenow = timenow; 01679 sb_threads[i].ilast = left; 01680 left = left - dec; 01681 if (left >0){ 01682 sb_threads[i].ifirst = left; 01683 } 01684 else 01685 sb_threads[i].ifirst = 0; 01686 sb_threads[i].do_effector = do_effector; 01687 sb_threads[i].do_deflector = 0;// not used here 01688 sb_threads[i].fieldfactor = 0.0f;// not used here 01689 sb_threads[i].windfactor = 0.0f;// not used here 01690 sb_threads[i].nr= i; 01691 sb_threads[i].tot= totthread; 01692 } 01693 if(totthread > 1) { 01694 BLI_init_threads(&threads, exec_scan_for_ext_spring_forces, totthread); 01695 01696 for(i=0; i<totthread; i++) 01697 BLI_insert_thread(&threads, &sb_threads[i]); 01698 01699 BLI_end_threads(&threads); 01700 } 01701 else 01702 exec_scan_for_ext_spring_forces(&sb_threads[0]); 01703 /* clean up */ 01704 MEM_freeN(sb_threads); 01705 01706 pdEndEffectors(&do_effector); 01707 } 01708 01709 01710 /* --- the spring external section*/ 01711 01712 static int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb,float*cc) 01713 { 01714 float mindist,cp; 01715 int winner =1; 01716 mindist = ABS(dot_v3v3(pos,a)); 01717 01718 cp = ABS(dot_v3v3(pos,b)); 01719 if ( mindist < cp ){ 01720 mindist = cp; 01721 winner =2; 01722 } 01723 01724 cp = ABS(dot_v3v3(pos,c)); 01725 if (mindist < cp ){ 01726 mindist = cp; 01727 winner =3; 01728 } 01729 switch (winner){ 01730 case 1: VECCOPY(w,ca); break; 01731 case 2: VECCOPY(w,cb); break; 01732 case 3: VECCOPY(w,cc); 01733 } 01734 return(winner); 01735 } 01736 01737 01738 01739 static int sb_detect_vertex_collisionCached(float opco[3], float facenormal[3], float *damp, 01740 float force[3], unsigned int UNUSED(par_layer), struct Object *vertexowner, 01741 float time,float vel[3], float *intrusion) 01742 { 01743 Object *ob= NULL; 01744 GHash *hash; 01745 GHashIterator *ihash; 01746 float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3],d_nvect[3], dv1[3],ve[3],avel[3]={0.0,0.0,0.0}, 01747 vv1[3], vv2[3], vv3[3], vv4[3], coledge[3]={0.0f, 0.0f, 0.0f}, mindistedge = 1000.0f, 01748 outerforceaccu[3],innerforceaccu[3], 01749 facedist,n_mag,force_mag_norm,minx,miny,minz,maxx,maxy,maxz, 01750 innerfacethickness = -0.5f, outerfacethickness = 0.2f, 01751 ee = 5.0f, ff = 0.1f, fa=1; 01752 int a, deflected=0, cavel=0,ci=0; 01753 /* init */ 01754 *intrusion = 0.0f; 01755 hash = vertexowner->soft->scratch->colliderhash; 01756 ihash = BLI_ghashIterator_new(hash); 01757 outerforceaccu[0]=outerforceaccu[1]=outerforceaccu[2]=0.0f; 01758 innerforceaccu[0]=innerforceaccu[1]=innerforceaccu[2]=0.0f; 01759 /* go */ 01760 while (!BLI_ghashIterator_isDone(ihash) ) { 01761 01762 ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); 01763 ob = BLI_ghashIterator_getKey (ihash); 01764 /* only with deflecting set */ 01765 if(ob->pd && ob->pd->deflect) { 01766 MFace *mface= NULL; 01767 MVert *mvert= NULL; 01768 MVert *mprevvert= NULL; 01769 ccdf_minmax *mima= NULL; 01770 01771 if(ccdm){ 01772 mface= ccdm->mface; 01773 mvert= ccdm->mvert; 01774 mprevvert= ccdm->mprevvert; 01775 mima= ccdm->mima; 01776 a = ccdm->totface; 01777 01778 minx =ccdm->bbmin[0]; 01779 miny =ccdm->bbmin[1]; 01780 minz =ccdm->bbmin[2]; 01781 01782 maxx =ccdm->bbmax[0]; 01783 maxy =ccdm->bbmax[1]; 01784 maxz =ccdm->bbmax[2]; 01785 01786 if ((opco[0] < minx) || 01787 (opco[1] < miny) || 01788 (opco[2] < minz) || 01789 (opco[0] > maxx) || 01790 (opco[1] > maxy) || 01791 (opco[2] > maxz) ) { 01792 /* outside the padded boundbox --> collision object is too far away */ 01793 BLI_ghashIterator_step(ihash); 01794 continue; 01795 } 01796 } 01797 else{ 01798 /*aye that should be cached*/ 01799 printf("missing cache error \n"); 01800 BLI_ghashIterator_step(ihash); 01801 continue; 01802 } 01803 01804 /* do object level stuff */ 01805 /* need to have user control for that since it depends on model scale */ 01806 innerfacethickness =-ob->pd->pdef_sbift; 01807 outerfacethickness =ob->pd->pdef_sboft; 01808 fa = (ff*outerfacethickness-outerfacethickness); 01809 fa *= fa; 01810 fa = 1.0f/fa; 01811 avel[0]=avel[1]=avel[2]=0.0f; 01812 /* use mesh*/ 01813 while (a--) { 01814 if ( 01815 (opco[0] < mima->minx) || 01816 (opco[0] > mima->maxx) || 01817 (opco[1] < mima->miny) || 01818 (opco[1] > mima->maxy) || 01819 (opco[2] < mima->minz) || 01820 (opco[2] > mima->maxz) 01821 ) { 01822 mface++; 01823 mima++; 01824 continue; 01825 } 01826 01827 if (mvert){ 01828 01829 VECCOPY(nv1,mvert[mface->v1].co); 01830 VECCOPY(nv2,mvert[mface->v2].co); 01831 VECCOPY(nv3,mvert[mface->v3].co); 01832 if (mface->v4){ 01833 VECCOPY(nv4,mvert[mface->v4].co); 01834 } 01835 01836 if (mprevvert){ 01837 /* grab the average speed of the collider vertices 01838 before we spoil nvX 01839 humm could be done once a SB steps but then we' need to store that too 01840 since the AABB reduced propabitlty to get here drasticallly 01841 it might be a nice tradeof CPU <--> memory 01842 */ 01843 VECSUB(vv1,nv1,mprevvert[mface->v1].co); 01844 VECSUB(vv2,nv2,mprevvert[mface->v2].co); 01845 VECSUB(vv3,nv3,mprevvert[mface->v3].co); 01846 if (mface->v4){ 01847 VECSUB(vv4,nv4,mprevvert[mface->v4].co); 01848 } 01849 01850 mul_v3_fl(nv1,time); 01851 Vec3PlusStVec(nv1,(1.0f-time),mprevvert[mface->v1].co); 01852 01853 mul_v3_fl(nv2,time); 01854 Vec3PlusStVec(nv2,(1.0f-time),mprevvert[mface->v2].co); 01855 01856 mul_v3_fl(nv3,time); 01857 Vec3PlusStVec(nv3,(1.0f-time),mprevvert[mface->v3].co); 01858 01859 if (mface->v4){ 01860 mul_v3_fl(nv4,time); 01861 Vec3PlusStVec(nv4,(1.0f-time),mprevvert[mface->v4].co); 01862 } 01863 } 01864 } 01865 01866 /* switch origin to be nv2*/ 01867 VECSUB(edge1, nv1, nv2); 01868 VECSUB(edge2, nv3, nv2); 01869 VECSUB(dv1,opco,nv2); /* abuse dv1 to have vertex in question at *origin* of triangle */ 01870 01871 cross_v3_v3v3(d_nvect, edge2, edge1); 01872 n_mag = normalize_v3(d_nvect); 01873 facedist = dot_v3v3(dv1,d_nvect); 01874 // so rules are 01875 // 01876 01877 if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){ 01878 if (isect_point_tri_prism_v3(opco, nv1, nv2, nv3) ){ 01879 force_mag_norm =(float)exp(-ee*facedist); 01880 if (facedist > outerfacethickness*ff) 01881 force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); 01882 *damp=ob->pd->pdef_sbdamp; 01883 if (facedist > 0.0f){ 01884 *damp *= (1.0f - facedist/outerfacethickness); 01885 Vec3PlusStVec(outerforceaccu,force_mag_norm,d_nvect); 01886 deflected = 3; 01887 01888 } 01889 else { 01890 Vec3PlusStVec(innerforceaccu,force_mag_norm,d_nvect); 01891 if (deflected < 2) deflected = 2; 01892 } 01893 if ((mprevvert) && (*damp > 0.0f)){ 01894 choose_winner(ve,opco,nv1,nv2,nv3,vv1,vv2,vv3); 01895 VECADD(avel,avel,ve); 01896 cavel ++; 01897 } 01898 *intrusion += facedist; 01899 ci++; 01900 } 01901 } 01902 if (mface->v4){ /* quad */ 01903 /* switch origin to be nv4 */ 01904 VECSUB(edge1, nv3, nv4); 01905 VECSUB(edge2, nv1, nv4); 01906 VECSUB(dv1,opco,nv4); /* abuse dv1 to have vertex in question at *origin* of triangle */ 01907 01908 cross_v3_v3v3(d_nvect, edge2, edge1); 01909 n_mag = normalize_v3(d_nvect); 01910 facedist = dot_v3v3(dv1,d_nvect); 01911 01912 if ((facedist > innerfacethickness) && (facedist < outerfacethickness)){ 01913 if (isect_point_tri_prism_v3(opco, nv1, nv3, nv4) ){ 01914 force_mag_norm =(float)exp(-ee*facedist); 01915 if (facedist > outerfacethickness*ff) 01916 force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); 01917 *damp=ob->pd->pdef_sbdamp; 01918 if (facedist > 0.0f){ 01919 *damp *= (1.0f - facedist/outerfacethickness); 01920 Vec3PlusStVec(outerforceaccu,force_mag_norm,d_nvect); 01921 deflected = 3; 01922 01923 } 01924 else { 01925 Vec3PlusStVec(innerforceaccu,force_mag_norm,d_nvect); 01926 if (deflected < 2) deflected = 2; 01927 } 01928 01929 if ((mprevvert) && (*damp > 0.0f)){ 01930 choose_winner(ve,opco,nv1,nv3,nv4,vv1,vv3,vv4); 01931 VECADD(avel,avel,ve); 01932 cavel ++; 01933 } 01934 *intrusion += facedist; 01935 ci++; 01936 } 01937 01938 } 01939 if ((deflected < 2)&& (G.rt != 444)) // we did not hit a face until now 01940 { // see if 'outer' hits an edge 01941 float dist; 01942 01943 closest_to_line_segment_v3(ve, opco, nv1, nv2); 01944 VECSUB(ve,opco,ve); 01945 dist = normalize_v3(ve); 01946 if ((dist < outerfacethickness)&&(dist < mindistedge )){ 01947 VECCOPY(coledge,ve); 01948 mindistedge = dist, 01949 deflected=1; 01950 } 01951 01952 closest_to_line_segment_v3(ve, opco, nv2, nv3); 01953 VECSUB(ve,opco,ve); 01954 dist = normalize_v3(ve); 01955 if ((dist < outerfacethickness)&&(dist < mindistedge )){ 01956 VECCOPY(coledge,ve); 01957 mindistedge = dist, 01958 deflected=1; 01959 } 01960 01961 closest_to_line_segment_v3(ve, opco, nv3, nv1); 01962 VECSUB(ve,opco,ve); 01963 dist = normalize_v3(ve); 01964 if ((dist < outerfacethickness)&&(dist < mindistedge )){ 01965 VECCOPY(coledge,ve); 01966 mindistedge = dist, 01967 deflected=1; 01968 } 01969 if (mface->v4){ /* quad */ 01970 closest_to_line_segment_v3(ve, opco, nv3, nv4); 01971 VECSUB(ve,opco,ve); 01972 dist = normalize_v3(ve); 01973 if ((dist < outerfacethickness)&&(dist < mindistedge )){ 01974 VECCOPY(coledge,ve); 01975 mindistedge = dist, 01976 deflected=1; 01977 } 01978 01979 closest_to_line_segment_v3(ve, opco, nv1, nv4); 01980 VECSUB(ve,opco,ve); 01981 dist = normalize_v3(ve); 01982 if ((dist < outerfacethickness)&&(dist < mindistedge )){ 01983 VECCOPY(coledge,ve); 01984 mindistedge = dist, 01985 deflected=1; 01986 } 01987 01988 } 01989 01990 01991 } 01992 } 01993 mface++; 01994 mima++; 01995 }/* while a */ 01996 } /* if(ob->pd && ob->pd->deflect) */ 01997 BLI_ghashIterator_step(ihash); 01998 } /* while () */ 01999 02000 if (deflected == 1){ // no face but 'outer' edge cylinder sees vert 02001 force_mag_norm =(float)exp(-ee*mindistedge); 02002 if (mindistedge > outerfacethickness*ff) 02003 force_mag_norm =(float)force_mag_norm*fa*(mindistedge - outerfacethickness)*(mindistedge - outerfacethickness); 02004 Vec3PlusStVec(force,force_mag_norm,coledge); 02005 *damp=ob->pd->pdef_sbdamp; 02006 if (mindistedge > 0.0f){ 02007 *damp *= (1.0f - mindistedge/outerfacethickness); 02008 } 02009 02010 } 02011 if (deflected == 2){ // face inner detected 02012 VECADD(force,force,innerforceaccu); 02013 } 02014 if (deflected == 3){ // face outer detected 02015 VECADD(force,force,outerforceaccu); 02016 } 02017 02018 BLI_ghashIterator_free(ihash); 02019 if (cavel) mul_v3_fl(avel,1.0f/(float)cavel); 02020 VECCOPY(vel,avel); 02021 if (ci) *intrusion /= ci; 02022 if (deflected){ 02023 normalize_v3_v3(facenormal, force); 02024 } 02025 return deflected; 02026 } 02027 02028 02029 /* sandbox to plug in various deflection algos */ 02030 static int sb_deflect_face(Object *ob,float *actpos,float *facenormal,float *force,float *cf,float time,float *vel,float *intrusion) 02031 { 02032 float s_actpos[3]; 02033 int deflected; 02034 VECCOPY(s_actpos,actpos); 02035 deflected= sb_detect_vertex_collisionCached(s_actpos, facenormal, cf, force , ob->lay, ob,time,vel,intrusion); 02036 //deflected= sb_detect_vertex_collisionCachedEx(s_actpos, facenormal, cf, force , ob->lay, ob,time,vel,intrusion); 02037 return(deflected); 02038 } 02039 02040 /* hiding this for now .. but the jacobian may pop up on other tasks .. so i'd like to keep it 02041 static void dfdx_spring(int ia, int ic, int op, float dir[3],float L,float len,float factor) 02042 { 02043 float m,delta_ij; 02044 int i ,j; 02045 if (L < len){ 02046 for(i=0;i<3;i++) 02047 for(j=0;j<3;j++){ 02048 delta_ij = (i==j ? (1.0f): (0.0f)); 02049 m=factor*(dir[i]*dir[j] + (1-L/len)*(delta_ij - dir[i]*dir[j])); 02050 nlMatrixAdd(ia+i,op+ic+j,m); 02051 } 02052 } 02053 else{ 02054 for(i=0;i<3;i++) 02055 for(j=0;j<3;j++){ 02056 m=factor*dir[i]*dir[j]; 02057 nlMatrixAdd(ia+i,op+ic+j,m); 02058 } 02059 } 02060 } 02061 02062 02063 static void dfdx_goal(int ia, int ic, int op, float factor) 02064 { 02065 int i; 02066 for(i=0;i<3;i++) nlMatrixAdd(ia+i,op+ic+i,factor); 02067 } 02068 02069 static void dfdv_goal(int ia, int ic,float factor) 02070 { 02071 int i; 02072 for(i=0;i<3;i++) nlMatrixAdd(ia+i,ic+i,factor); 02073 } 02074 */ 02075 static void sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float UNUSED(forcetime), int nl_flags) 02076 { 02077 SoftBody *sb= ob->soft; /* is supposed to be there */ 02078 BodyPoint *bp1,*bp2; 02079 02080 float dir[3],dvel[3]; 02081 float distance,forcefactor,kd,absvel,projvel,kw; 02082 #if 0 /* UNUSED */ 02083 int ia,ic; 02084 #endif 02085 /* prepare depending on which side of the spring we are on */ 02086 if (bpi == bs->v1){ 02087 bp1 = &sb->bpoint[bs->v1]; 02088 bp2 = &sb->bpoint[bs->v2]; 02089 #if 0 /* UNUSED */ 02090 ia =3*bs->v1; 02091 ic =3*bs->v2; 02092 #endif 02093 } 02094 else if (bpi == bs->v2){ 02095 bp1 = &sb->bpoint[bs->v2]; 02096 bp2 = &sb->bpoint[bs->v1]; 02097 #if 0 /* UNUSED */ 02098 ia =3*bs->v2; 02099 ic =3*bs->v1; 02100 #endif 02101 } 02102 else{ 02103 /* TODO make this debug option */ 02104 02105 printf("bodypoint <bpi> is not attached to spring <*bs> --> sb_spring_force()\n"); 02106 return; 02107 } 02108 02109 /* do bp1 <--> bp2 elastic */ 02110 sub_v3_v3v3(dir,bp1->pos,bp2->pos); 02111 distance = normalize_v3(dir); 02112 if (bs->len < distance) 02113 iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ 02114 else 02115 iks = 1.0f/(1.0f-sb->inpush)-1.0f ;/* inner spring constants function */ 02116 02117 if(bs->len > 0.0f) /* check for degenerated springs */ 02118 forcefactor = iks/bs->len; 02119 else 02120 forcefactor = iks; 02121 kw = (bp1->springweight+bp2->springweight)/2.0f; 02122 kw = kw * kw; 02123 kw = kw * kw; 02124 switch (bs->springtype){ 02125 case SB_EDGE: 02126 case SB_HANDLE: 02127 forcefactor *= kw; 02128 break; 02129 case SB_BEND: 02130 forcefactor *=sb->secondspring*kw; 02131 break; 02132 case SB_STIFFQUAD: 02133 forcefactor *=sb->shearstiff*sb->shearstiff* kw; 02134 break; 02135 default: 02136 break; 02137 } 02138 02139 02140 Vec3PlusStVec(bp1->force,(bs->len - distance)*forcefactor,dir); 02141 02142 /* do bp1 <--> bp2 viscous */ 02143 sub_v3_v3v3(dvel,bp1->vec,bp2->vec); 02144 kd = sb->infrict * sb_fric_force_scale(ob); 02145 absvel = normalize_v3(dvel); 02146 projvel = dot_v3v3(dir,dvel); 02147 kd *= absvel * projvel; 02148 Vec3PlusStVec(bp1->force,-kd,dir); 02149 02150 /* do jacobian stuff if needed */ 02151 if(nl_flags & NLF_BUILD){ 02152 //int op =3*sb->totpoint; 02153 //float mvel = -forcetime*kd; 02154 //float mpos = -forcetime*forcefactor; 02155 /* depending on my pos */ 02156 // dfdx_spring(ia,ia,op,dir,bs->len,distance,-mpos); 02157 /* depending on my vel */ 02158 // dfdv_goal(ia,ia,mvel); // well that ignores geometie 02159 if(bp2->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */ 02160 /* depending on other pos */ 02161 // dfdx_spring(ia,ic,op,dir,bs->len,distance,mpos); 02162 /* depending on other vel */ 02163 // dfdv_goal(ia,ia,-mvel); // well that ignores geometie 02164 } 02165 } 02166 } 02167 02168 02169 /* since this is definitely the most CPU consuming task here .. try to spread it */ 02170 /* core function _softbody_calc_forces_slice_in_a_thread */ 02171 /* result is int to be able to flag user break */ 02172 static int _softbody_calc_forces_slice_in_a_thread(Scene *scene, Object *ob, float forcetime, float timenow,int ifirst,int ilast,int *UNUSED(ptr_to_break_func(void)),ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor) 02173 { 02174 float iks; 02175 int bb,do_selfcollision,do_springcollision,do_aero; 02176 int number_of_points_here = ilast - ifirst; 02177 SoftBody *sb= ob->soft; /* is supposed to be there */ 02178 BodyPoint *bp; 02179 02180 /* intitialize */ 02181 if (sb) { 02182 /* check conditions for various options */ 02183 /* +++ could be done on object level to squeeze out the last bits of it */ 02184 do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); 02185 do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); 02186 do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); 02187 /* --- could be done on object level to squeeze out the last bits of it */ 02188 } 02189 else { 02190 printf("Error expected a SB here \n"); 02191 return (999); 02192 } 02193 02194 /* debugerin */ 02195 if (sb->totpoint < ifirst) { 02196 printf("Aye 998"); 02197 return (998); 02198 } 02199 /* debugerin */ 02200 02201 02202 bp = &sb->bpoint[ifirst]; 02203 for(bb=number_of_points_here; bb>0; bb--, bp++) { 02204 /* clear forces accumulator */ 02205 bp->force[0]= bp->force[1]= bp->force[2]= 0.0; 02206 /* naive ball self collision */ 02207 /* needs to be done if goal snaps or not */ 02208 if(do_selfcollision){ 02209 int attached; 02210 BodyPoint *obp; 02211 BodySpring *bs; 02212 int c,b; 02213 float velcenter[3],dvel[3],def[3]; 02214 float distance; 02215 float compare; 02216 float bstune = sb->ballstiff; 02217 02218 for(c=sb->totpoint, obp= sb->bpoint; c>=ifirst+bb; c--, obp++) { 02219 compare = (obp->colball + bp->colball); 02220 sub_v3_v3v3(def, bp->pos, obp->pos); 02221 /* rather check the AABBoxes before ever calulating the real distance */ 02222 /* mathematically it is completly nuts, but performace is pretty much (3) times faster */ 02223 if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; 02224 distance = normalize_v3(def); 02225 if (distance < compare ){ 02226 /* exclude body points attached with a spring */ 02227 attached = 0; 02228 for(b=obp->nofsprings;b>0;b--){ 02229 bs = sb->bspring + obp->springs[b-1]; 02230 if (( ilast-bb == bs->v2) || ( ilast-bb == bs->v1)){ 02231 attached=1; 02232 continue;} 02233 } 02234 if (!attached){ 02235 float f = bstune/(distance) + bstune/(compare*compare)*distance - 2.0f*bstune/compare ; 02236 02237 mid_v3_v3v3(velcenter, bp->vec, obp->vec); 02238 sub_v3_v3v3(dvel,velcenter,bp->vec); 02239 mul_v3_fl(dvel,_final_mass(ob,bp)); 02240 02241 Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def); 02242 Vec3PlusStVec(bp->force,sb->balldamp,dvel); 02243 02244 /* exploit force(a,b) == -force(b,a) part2/2 */ 02245 sub_v3_v3v3(dvel,velcenter,obp->vec); 02246 mul_v3_fl(dvel,_final_mass(ob,bp)); 02247 02248 Vec3PlusStVec(obp->force,sb->balldamp,dvel); 02249 Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def); 02250 } 02251 } 02252 } 02253 } 02254 /* naive ball self collision done */ 02255 02256 if(_final_goal(ob,bp) < SOFTGOALSNAP){ /* ommit this bp when it snaps */ 02257 float auxvect[3]; 02258 float velgoal[3]; 02259 02260 /* do goal stuff */ 02261 if(ob->softflag & OB_SB_GOAL) { 02262 /* true elastic goal */ 02263 float ks,kd; 02264 sub_v3_v3v3(auxvect,bp->pos,bp->origT); 02265 ks = 1.0f/(1.0f- _final_goal(ob,bp)*sb->goalspring)-1.0f ; 02266 bp->force[0]+= -ks*(auxvect[0]); 02267 bp->force[1]+= -ks*(auxvect[1]); 02268 bp->force[2]+= -ks*(auxvect[2]); 02269 02270 /* calulate damping forces generated by goals*/ 02271 sub_v3_v3v3(velgoal,bp->origS, bp->origE); 02272 kd = sb->goalfrict * sb_fric_force_scale(ob) ; 02273 add_v3_v3v3(auxvect,velgoal,bp->vec); 02274 02275 if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */ 02276 bp->force[0]-= kd * (auxvect[0]); 02277 bp->force[1]-= kd * (auxvect[1]); 02278 bp->force[2]-= kd * (auxvect[2]); 02279 } 02280 else { 02281 bp->force[0]-= kd * (velgoal[0] - bp->vec[0]); 02282 bp->force[1]-= kd * (velgoal[1] - bp->vec[1]); 02283 bp->force[2]-= kd * (velgoal[2] - bp->vec[2]); 02284 } 02285 } 02286 /* done goal stuff */ 02287 02288 /* gravitation */ 02289 if (sb && scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY){ 02290 float gravity[3]; 02291 VECCOPY(gravity, scene->physics_settings.gravity); 02292 mul_v3_fl(gravity, sb_grav_force_scale(ob)*_final_mass(ob,bp)*sb->effector_weights->global_gravity); /* individual mass of node here */ 02293 add_v3_v3(bp->force, gravity); 02294 } 02295 02296 /* particle field & vortex */ 02297 if(do_effector) { 02298 EffectedPoint epoint; 02299 float kd; 02300 float force[3]= {0.0f, 0.0f, 0.0f}; 02301 float speed[3]= {0.0f, 0.0f, 0.0f}; 02302 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ 02303 pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint); 02304 pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); 02305 02306 /* apply forcefield*/ 02307 mul_v3_fl(force,fieldfactor* eval_sb_fric_force_scale); 02308 VECADD(bp->force, bp->force, force); 02309 02310 /* BP friction in moving media */ 02311 kd= sb->mediafrict* eval_sb_fric_force_scale; 02312 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale); 02313 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale); 02314 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale); 02315 /* now we'll have nice centrifugal effect for vortex */ 02316 02317 } 02318 else { 02319 /* BP friction in media (not) moving*/ 02320 float kd = sb->mediafrict* sb_fric_force_scale(ob); 02321 /* assume it to be proportional to actual velocity */ 02322 bp->force[0]-= bp->vec[0]*kd; 02323 bp->force[1]-= bp->vec[1]*kd; 02324 bp->force[2]-= bp->vec[2]*kd; 02325 /* friction in media done */ 02326 } 02327 /* +++cached collision targets */ 02328 bp->choke = 0.0f; 02329 bp->choke2 = 0.0f; 02330 bp->loc_flag &= ~SBF_DOFUZZY; 02331 if(do_deflector && !(bp->loc_flag & SBF_OUTOFCOLLISION) ) { 02332 float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion; 02333 float kd = 1.0f; 02334 02335 if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){ 02336 if (intrusion < 0.0f){ 02337 sb->scratch->flag |= SBF_DOFUZZY; 02338 bp->loc_flag |= SBF_DOFUZZY; 02339 bp->choke = sb->choke*0.01f; 02340 } 02341 02342 VECSUB(cfforce,bp->vec,vel); 02343 Vec3PlusStVec(bp->force,-cf*50.0f,cfforce); 02344 02345 Vec3PlusStVec(bp->force,kd,defforce); 02346 } 02347 02348 } 02349 /* ---cached collision targets */ 02350 02351 /* +++springs */ 02352 iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ 02353 if(ob->softflag & OB_SB_EDGES) { 02354 if (sb->bspring){ /* spring list exists at all ? */ 02355 int b; 02356 BodySpring *bs; 02357 for(b=bp->nofsprings;b>0;b--){ 02358 bs = sb->bspring + bp->springs[b-1]; 02359 if (do_springcollision || do_aero){ 02360 add_v3_v3(bp->force, bs->ext_force); 02361 if (bs->flag & BSF_INTERSECT) 02362 bp->choke = bs->cf; 02363 02364 } 02365 // sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags) 02366 sb_spring_force(ob,ilast-bb,bs,iks,forcetime,0); 02367 }/* loop springs */ 02368 }/* existing spring list */ 02369 }/*any edges*/ 02370 /* ---springs */ 02371 }/*omit on snap */ 02372 }/*loop all bp's*/ 02373 return 0; /*done fine*/ 02374 } 02375 02376 static void *exec_softbody_calc_forces(void *data) 02377 { 02378 SB_thread_context *pctx = (SB_thread_context*)data; 02379 _softbody_calc_forces_slice_in_a_thread(pctx->scene, pctx->ob, pctx->forcetime, pctx->timenow, pctx->ifirst, pctx->ilast, NULL, pctx->do_effector,pctx->do_deflector,pctx->fieldfactor,pctx->windfactor); 02380 return NULL; 02381 } 02382 02383 static void sb_cf_threads_run(Scene *scene, Object *ob, float forcetime, float timenow,int totpoint,int *UNUSED(ptr_to_break_func(void)),struct ListBase *do_effector,int do_deflector,float fieldfactor, float windfactor) 02384 { 02385 ListBase threads; 02386 SB_thread_context *sb_threads; 02387 int i, totthread,left,dec; 02388 int lowpoints =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ 02389 02390 /* figure the number of threads while preventing pretty pointless threading overhead */ 02391 if(scene->r.mode & R_FIXED_THREADS) 02392 totthread= scene->r.threads; 02393 else 02394 totthread= BLI_system_thread_count(); 02395 /* what if we got zillions of CPUs running but less to spread*/ 02396 while ((totpoint/totthread < lowpoints) && (totthread > 1)) { 02397 totthread--; 02398 } 02399 02400 /* printf("sb_cf_threads_run spawning %d threads \n",totthread); */ 02401 02402 sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBThread"); 02403 memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); 02404 left = totpoint; 02405 dec = totpoint/totthread +1; 02406 for(i=0; i<totthread; i++) { 02407 sb_threads[i].scene = scene; 02408 sb_threads[i].ob = ob; 02409 sb_threads[i].forcetime = forcetime; 02410 sb_threads[i].timenow = timenow; 02411 sb_threads[i].ilast = left; 02412 left = left - dec; 02413 if (left >0){ 02414 sb_threads[i].ifirst = left; 02415 } 02416 else 02417 sb_threads[i].ifirst = 0; 02418 sb_threads[i].do_effector = do_effector; 02419 sb_threads[i].do_deflector = do_deflector; 02420 sb_threads[i].fieldfactor = fieldfactor; 02421 sb_threads[i].windfactor = windfactor; 02422 sb_threads[i].nr= i; 02423 sb_threads[i].tot= totthread; 02424 } 02425 02426 02427 if(totthread > 1) { 02428 BLI_init_threads(&threads, exec_softbody_calc_forces, totthread); 02429 02430 for(i=0; i<totthread; i++) 02431 BLI_insert_thread(&threads, &sb_threads[i]); 02432 02433 BLI_end_threads(&threads); 02434 } 02435 else 02436 exec_softbody_calc_forces(&sb_threads[0]); 02437 /* clean up */ 02438 MEM_freeN(sb_threads); 02439 } 02440 02441 static void softbody_calc_forcesEx(Scene *scene, Object *ob, float forcetime, float timenow, int UNUSED(nl_flags)) 02442 { 02443 /* rule we never alter free variables :bp->vec bp->pos in here ! 02444 * this will ruin adaptive stepsize AKA heun! (BM) 02445 */ 02446 SoftBody *sb= ob->soft; /* is supposed to be there */ 02447 /*BodyPoint *bproot;*/ /* UNUSED */ 02448 ListBase *do_effector = NULL; 02449 /* float gravity; */ /* UNUSED */ 02450 /* float iks; */ 02451 float fieldfactor = -1.0f, windfactor = 0.25; 02452 int do_deflector /*,do_selfcollision*/ ,do_springcollision,do_aero; 02453 02454 /* gravity = sb->grav * sb_grav_force_scale(ob); */ /* UNUSED */ 02455 02456 /* check conditions for various options */ 02457 do_deflector= query_external_colliders(scene, ob); 02458 /* do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); */ /* UNUSED */ 02459 do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); 02460 do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); 02461 02462 /* iks = 1.0f/(1.0f-sb->inspring)-1.0f; */ /* inner spring constants function */ /* UNUSED */ 02463 /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */ 02464 02465 if (do_springcollision || do_aero) 02466 sb_sfesf_threads_run(scene, ob, timenow,sb->totspring,NULL); 02467 02468 /* after spring scan because it uses Effoctors too */ 02469 do_effector= pdInitEffectors(scene, ob, NULL, sb->effector_weights); 02470 02471 if (do_deflector) { 02472 float defforce[3]; 02473 do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow); 02474 } 02475 02476 sb_cf_threads_run(scene, ob, forcetime, timenow, sb->totpoint, NULL, do_effector, do_deflector, fieldfactor, windfactor); 02477 02478 /* finally add forces caused by face collision */ 02479 if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow); 02480 02481 /* finish matrix and solve */ 02482 pdEndEffectors(&do_effector); 02483 } 02484 02485 02486 02487 02488 static void softbody_calc_forces(Scene *scene, Object *ob, float forcetime, float timenow, int nl_flags) 02489 { 02490 /* redirection to the new threaded Version */ 02491 if (!(G.rt & 0x10)){ // 16 02492 softbody_calc_forcesEx(scene, ob, forcetime, timenow, nl_flags); 02493 return; 02494 } 02495 else{ 02496 /* so the following will die */ 02497 /* |||||||||||||||||||||||||| */ 02498 /* VVVVVVVVVVVVVVVVVVVVVVVVVV */ 02499 /*backward compatibility note: 02500 fixing bug [17428] which forces adaptive step size to tiny steps 02501 in some situations 02502 .. keeping G.rt==17 0x11 option for old files 'needing' the bug*/ 02503 02504 /* rule we never alter free variables :bp->vec bp->pos in here ! 02505 * this will ruin adaptive stepsize AKA heun! (BM) 02506 */ 02507 SoftBody *sb= ob->soft; /* is supposed to be there */ 02508 BodyPoint *bp; 02509 /* BodyPoint *bproot; */ /* UNUSED */ 02510 BodySpring *bs; 02511 ListBase *do_effector = NULL; 02512 float iks, ks, kd, gravity[3] = {0.0f,0.0f,0.0f}; 02513 float fieldfactor = -1.0f, windfactor = 0.25f; 02514 float tune = sb->ballstiff; 02515 int a, b, do_deflector,do_selfcollision,do_springcollision,do_aero; 02516 02517 02518 /* jacobian 02519 NLboolean success; 02520 02521 if(nl_flags){ 02522 nlBegin(NL_SYSTEM); 02523 nlBegin(NL_MATRIX); 02524 } 02525 */ 02526 02527 02528 if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY){ 02529 VECCOPY(gravity, scene->physics_settings.gravity); 02530 mul_v3_fl(gravity, sb_grav_force_scale(ob)*sb->effector_weights->global_gravity); 02531 } 02532 02533 /* check conditions for various options */ 02534 do_deflector= query_external_colliders(scene, ob); 02535 do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); 02536 do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); 02537 do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); 02538 02539 iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ 02540 /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */ 02541 02542 if (do_springcollision || do_aero) scan_for_ext_spring_forces(scene, ob, timenow); 02543 /* after spring scan because it uses Effoctors too */ 02544 do_effector= pdInitEffectors(scene, ob, NULL, ob->soft->effector_weights); 02545 02546 if (do_deflector) { 02547 float defforce[3]; 02548 do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow); 02549 } 02550 02551 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02552 /* clear forces accumulator */ 02553 bp->force[0]= bp->force[1]= bp->force[2]= 0.0; 02554 if(nl_flags & NLF_BUILD){ 02555 //int ia =3*(sb->totpoint-a); 02556 //int op =3*sb->totpoint; 02557 /* dF/dV = v */ 02558 /* jacobioan 02559 nlMatrixAdd(op+ia,ia,-forcetime); 02560 nlMatrixAdd(op+ia+1,ia+1,-forcetime); 02561 nlMatrixAdd(op+ia+2,ia+2,-forcetime); 02562 02563 nlMatrixAdd(ia,ia,1); 02564 nlMatrixAdd(ia+1,ia+1,1); 02565 nlMatrixAdd(ia+2,ia+2,1); 02566 02567 nlMatrixAdd(op+ia,op+ia,1); 02568 nlMatrixAdd(op+ia+1,op+ia+1,1); 02569 nlMatrixAdd(op+ia+2,op+ia+2,1); 02570 */ 02571 02572 02573 } 02574 02575 /* naive ball self collision */ 02576 /* needs to be done if goal snaps or not */ 02577 if(do_selfcollision){ 02578 int attached; 02579 BodyPoint *obp; 02580 int c,b; 02581 float velcenter[3],dvel[3],def[3]; 02582 float distance; 02583 float compare; 02584 02585 for(c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) { 02586 02587 //if ((bp->octantflag & obp->octantflag) == 0) continue; 02588 02589 compare = (obp->colball + bp->colball); 02590 sub_v3_v3v3(def, bp->pos, obp->pos); 02591 02592 /* rather check the AABBoxes before ever calulating the real distance */ 02593 /* mathematically it is completly nuts, but performace is pretty much (3) times faster */ 02594 if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; 02595 02596 distance = normalize_v3(def); 02597 if (distance < compare ){ 02598 /* exclude body points attached with a spring */ 02599 attached = 0; 02600 for(b=obp->nofsprings;b>0;b--){ 02601 bs = sb->bspring + obp->springs[b-1]; 02602 if (( sb->totpoint-a == bs->v2) || ( sb->totpoint-a == bs->v1)){ 02603 attached=1; 02604 continue;} 02605 } 02606 if (!attached){ 02607 float f = tune/(distance) + tune/(compare*compare)*distance - 2.0f*tune/compare ; 02608 02609 mid_v3_v3v3(velcenter, bp->vec, obp->vec); 02610 sub_v3_v3v3(dvel,velcenter,bp->vec); 02611 mul_v3_fl(dvel,_final_mass(ob,bp)); 02612 02613 Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def); 02614 Vec3PlusStVec(bp->force,sb->balldamp,dvel); 02615 02616 if(nl_flags & NLF_BUILD){ 02617 //int ia =3*(sb->totpoint-a); 02618 //int ic =3*(sb->totpoint-c); 02619 //int op =3*sb->totpoint; 02620 //float mvel = forcetime*sb->nodemass*sb->balldamp; 02621 //float mpos = forcetime*tune*(1.0f-sb->balldamp); 02622 /*some quick and dirty entries to the jacobian*/ 02623 //dfdx_goal(ia,ia,op,mpos); 02624 //dfdv_goal(ia,ia,mvel); 02625 /* exploit force(a,b) == -force(b,a) part1/2 */ 02626 //dfdx_goal(ic,ic,op,mpos); 02627 //dfdv_goal(ic,ic,mvel); 02628 02629 02630 /*TODO sit down an X-out the true jacobian entries*/ 02631 /*well does not make to much sense because the eigenvalues 02632 of the jacobian go negative; and negative eigenvalues 02633 on a complex iterative system z(n+1)=A * z(n) 02634 give imaginary roots in the charcateristic polynom 02635 --> solutions that to z(t)=u(t)* exp ( i omega t) --> oscilations we don't want here 02636 where u(t) is a unknown amplitude function (worst case rising fast) 02637 */ 02638 } 02639 02640 /* exploit force(a,b) == -force(b,a) part2/2 */ 02641 sub_v3_v3v3(dvel,velcenter,obp->vec); 02642 mul_v3_fl(dvel,(_final_mass(ob,bp)+_final_mass(ob,obp))/2.0f); 02643 02644 Vec3PlusStVec(obp->force,sb->balldamp,dvel); 02645 Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def); 02646 02647 02648 } 02649 } 02650 } 02651 } 02652 /* naive ball self collision done */ 02653 02654 if(_final_goal(ob,bp) < SOFTGOALSNAP){ /* ommit this bp when it snaps */ 02655 float auxvect[3]; 02656 float velgoal[3]; 02657 02658 /* do goal stuff */ 02659 if(ob->softflag & OB_SB_GOAL) { 02660 /* true elastic goal */ 02661 sub_v3_v3v3(auxvect,bp->pos,bp->origT); 02662 ks = 1.0f/(1.0f- _final_goal(ob,bp)*sb->goalspring)-1.0f ; 02663 bp->force[0]+= -ks*(auxvect[0]); 02664 bp->force[1]+= -ks*(auxvect[1]); 02665 bp->force[2]+= -ks*(auxvect[2]); 02666 02667 if(nl_flags & NLF_BUILD){ 02668 //int ia =3*(sb->totpoint-a); 02669 //int op =3*(sb->totpoint); 02670 /* depending on my pos */ 02671 //dfdx_goal(ia,ia,op,ks*forcetime); 02672 } 02673 02674 02675 /* calulate damping forces generated by goals*/ 02676 sub_v3_v3v3(velgoal,bp->origS, bp->origE); 02677 kd = sb->goalfrict * sb_fric_force_scale(ob) ; 02678 add_v3_v3v3(auxvect,velgoal,bp->vec); 02679 02680 if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */ 02681 bp->force[0]-= kd * (auxvect[0]); 02682 bp->force[1]-= kd * (auxvect[1]); 02683 bp->force[2]-= kd * (auxvect[2]); 02684 if(nl_flags & NLF_BUILD){ 02685 //int ia =3*(sb->totpoint-a); 02686 normalize_v3(auxvect); 02687 /* depending on my vel */ 02688 //dfdv_goal(ia,ia,kd*forcetime); 02689 } 02690 02691 } 02692 else { 02693 bp->force[0]-= kd * (velgoal[0] - bp->vec[0]); 02694 bp->force[1]-= kd * (velgoal[1] - bp->vec[1]); 02695 bp->force[2]-= kd * (velgoal[2] - bp->vec[2]); 02696 } 02697 } 02698 /* done goal stuff */ 02699 02700 02701 /* gravitation */ 02702 VECADDFAC(bp->force, bp->force, gravity, _final_mass(ob,bp)); /* individual mass of node here */ 02703 02704 02705 /* particle field & vortex */ 02706 if(do_effector) { 02707 EffectedPoint epoint; 02708 float force[3]= {0.0f, 0.0f, 0.0f}; 02709 float speed[3]= {0.0f, 0.0f, 0.0f}; 02710 float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ 02711 pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint); 02712 pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); 02713 02714 /* apply forcefield*/ 02715 mul_v3_fl(force,fieldfactor* eval_sb_fric_force_scale); 02716 VECADD(bp->force, bp->force, force); 02717 02718 /* BP friction in moving media */ 02719 kd= sb->mediafrict* eval_sb_fric_force_scale; 02720 bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale); 02721 bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale); 02722 bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale); 02723 /* now we'll have nice centrifugal effect for vortex */ 02724 02725 } 02726 else { 02727 /* BP friction in media (not) moving*/ 02728 kd= sb->mediafrict* sb_fric_force_scale(ob); 02729 /* assume it to be proportional to actual velocity */ 02730 bp->force[0]-= bp->vec[0]*kd; 02731 bp->force[1]-= bp->vec[1]*kd; 02732 bp->force[2]-= bp->vec[2]*kd; 02733 /* friction in media done */ 02734 if(nl_flags & NLF_BUILD){ 02735 //int ia =3*(sb->totpoint-a); 02736 /* da/dv = */ 02737 02738 // nlMatrixAdd(ia,ia,forcetime*kd); 02739 // nlMatrixAdd(ia+1,ia+1,forcetime*kd); 02740 // nlMatrixAdd(ia+2,ia+2,forcetime*kd); 02741 } 02742 02743 } 02744 /* +++cached collision targets */ 02745 bp->choke = 0.0f; 02746 bp->choke2 = 0.0f; 02747 bp->loc_flag &= ~SBF_DOFUZZY; 02748 if(do_deflector) { 02749 float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion; 02750 kd = 1.0f; 02751 02752 if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){ 02753 if ((!nl_flags)&&(intrusion < 0.0f)){ 02754 if(G.rt & 0x01){ // 17 we did check for bit 0x10 before 02755 /*fixing bug [17428] this forces adaptive step size to tiny steps 02756 in some situations .. keeping G.rt==17 option for old files 'needing' the bug 02757 */ 02758 /*bjornmose: uugh.. what an evil hack 02759 violation of the 'don't touch bp->pos in here' rule 02760 but works nice, like this--> 02761 we predict the solution beeing out of the collider 02762 in heun step No1 and leave the heun step No2 adapt to it 02763 so we kind of introduced a implicit solver for this case 02764 */ 02765 Vec3PlusStVec(bp->pos,-intrusion,facenormal); 02766 } 02767 else{ 02768 02769 VECSUB(cfforce,bp->vec,vel); 02770 Vec3PlusStVec(bp->force,-cf*50.0f,cfforce); 02771 } 02772 02773 02774 sb->scratch->flag |= SBF_DOFUZZY; 02775 bp->loc_flag |= SBF_DOFUZZY; 02776 bp->choke = sb->choke*0.01f; 02777 } 02778 else{ 02779 VECSUB(cfforce,bp->vec,vel); 02780 Vec3PlusStVec(bp->force,-cf*50.0f,cfforce); 02781 } 02782 Vec3PlusStVec(bp->force,kd,defforce); 02783 if (nl_flags & NLF_BUILD){ 02784 // int ia =3*(sb->totpoint-a); 02785 // int op =3*sb->totpoint; 02786 //dfdx_goal(ia,ia,op,mpos); // don't do unless you know 02787 //dfdv_goal(ia,ia,-cf); 02788 02789 } 02790 02791 } 02792 02793 } 02794 /* ---cached collision targets */ 02795 02796 /* +++springs */ 02797 if(ob->softflag & OB_SB_EDGES) { 02798 if (sb->bspring){ /* spring list exists at all ? */ 02799 for(b=bp->nofsprings;b>0;b--){ 02800 bs = sb->bspring + bp->springs[b-1]; 02801 if (do_springcollision || do_aero){ 02802 add_v3_v3(bp->force, bs->ext_force); 02803 if (bs->flag & BSF_INTERSECT) 02804 bp->choke = bs->cf; 02805 02806 } 02807 // sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags) 02808 // rather remove nl_falgs from code .. will make things a lot cleaner 02809 sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,0); 02810 }/* loop springs */ 02811 }/* existing spring list */ 02812 }/*any edges*/ 02813 /* ---springs */ 02814 }/*omit on snap */ 02815 }/*loop all bp's*/ 02816 02817 02818 /* finally add forces caused by face collision */ 02819 if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow); 02820 02821 /* finish matrix and solve */ 02822 #if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be usefull .. so keep that BM 02823 if(nl_flags & NLF_SOLVE){ 02824 //double sct,sst=PIL_check_seconds_timer(); 02825 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02826 int iv =3*(sb->totpoint-a); 02827 int ip =3*(2*sb->totpoint-a); 02828 int n; 02829 for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);} 02830 for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);} 02831 } 02832 nlEnd(NL_MATRIX); 02833 nlEnd(NL_SYSTEM); 02834 02835 if ((G.rt == 32) && (nl_flags & NLF_BUILD)) 02836 { 02837 printf("####MEE#####\n"); 02838 nlPrintMatrix(); 02839 } 02840 02841 success= nlSolveAdvanced(NULL, 1); 02842 02843 // nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */ 02844 if(success){ 02845 float f; 02846 int index =0; 02847 /* for debug purpose .. anyhow cropping B vector looks like working */ 02848 if (G.rt ==32) 02849 for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02850 f=nlGetVariable(0,index); 02851 printf("(%f ",f);index++; 02852 f=nlGetVariable(0,index); 02853 printf("%f ",f);index++; 02854 f=nlGetVariable(0,index); 02855 printf("%f)",f);index++; 02856 } 02857 02858 index =0; 02859 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02860 f=nlGetVariable(0,index); 02861 bp->impdv[0] = f; index++; 02862 f=nlGetVariable(0,index); 02863 bp->impdv[1] = f; index++; 02864 f=nlGetVariable(0,index); 02865 bp->impdv[2] = f; index++; 02866 } 02867 /* 02868 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02869 f=nlGetVariable(0,index); 02870 bp->impdx[0] = f; index++; 02871 f=nlGetVariable(0,index); 02872 bp->impdx[1] = f; index++; 02873 f=nlGetVariable(0,index); 02874 bp->impdx[2] = f; index++; 02875 } 02876 */ 02877 } 02878 else{ 02879 printf("Matrix inversion failed \n"); 02880 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02881 VECCOPY(bp->impdv,bp->force); 02882 } 02883 02884 } 02885 02886 //sct=PIL_check_seconds_timer(); 02887 //if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name); 02888 } 02889 /* cleanup */ 02890 #endif 02891 pdEndEffectors(&do_effector); 02892 } 02893 } 02894 02895 static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags) 02896 { 02897 /* time evolution */ 02898 /* actually does an explicit euler step mode == 0 */ 02899 /* or heun ~ 2nd order runge-kutta steps, mode 1,2 */ 02900 SoftBody *sb= ob->soft; /* is supposed to be there */ 02901 BodyPoint *bp; 02902 float dx[3]={0},dv[3],aabbmin[3],aabbmax[3],cm[3]={0.0f,0.0f,0.0f}; 02903 float timeovermass/*,freezeloc=0.00001f,freezeforce=0.00000000001f*/; 02904 float maxerrpos= 0.0f,maxerrvel = 0.0f; 02905 int a,fuzzy=0; 02906 02907 forcetime *= sb_time_scale(ob); 02908 02909 aabbmin[0]=aabbmin[1]=aabbmin[2] = 1e20f; 02910 aabbmax[0]=aabbmax[1]=aabbmax[2] = -1e20f; 02911 02912 /* old one with homogenous masses */ 02913 /* claim a minimum mass for vertex */ 02914 /* 02915 if (sb->nodemass > 0.009999f) timeovermass = forcetime/sb->nodemass; 02916 else timeovermass = forcetime/0.009999f; 02917 */ 02918 02919 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 02920 /* now we have individual masses */ 02921 /* claim a minimum mass for vertex */ 02922 if (_final_mass(ob,bp) > 0.009999f) timeovermass = forcetime/_final_mass(ob,bp); 02923 else timeovermass = forcetime/0.009999f; 02924 02925 02926 if(_final_goal(ob,bp) < SOFTGOALSNAP){ 02927 /* this makes t~ = t */ 02928 if(mid_flags & MID_PRESERVE) VECCOPY(dx,bp->vec); 02929 02930 /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces + more forces*/ 02931 /* the ( ... )' operator denotes derivate respective time */ 02932 /* the euler step for velocity then becomes */ 02933 /* v(t + dt) = v(t) + a(t) * dt */ 02934 mul_v3_fl(bp->force,timeovermass);/* individual mass of node here */ 02935 /* some nasty if's to have heun in here too */ 02936 VECCOPY(dv,bp->force); 02937 02938 if (mode == 1){ 02939 VECCOPY(bp->prevvec, bp->vec); 02940 VECCOPY(bp->prevdv, dv); 02941 } 02942 02943 if (mode ==2){ 02944 /* be optimistic and execute step */ 02945 bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]); 02946 bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]); 02947 bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]); 02948 /* compare euler to heun to estimate error for step sizing */ 02949 maxerrvel = MAX2(maxerrvel,ABS(dv[0] - bp->prevdv[0])); 02950 maxerrvel = MAX2(maxerrvel,ABS(dv[1] - bp->prevdv[1])); 02951 maxerrvel = MAX2(maxerrvel,ABS(dv[2] - bp->prevdv[2])); 02952 } 02953 else {VECADD(bp->vec, bp->vec, bp->force);} 02954 02955 /* this makes t~ = t+dt */ 02956 if(!(mid_flags & MID_PRESERVE)) VECCOPY(dx,bp->vec); 02957 02958 /* so here is (x)'= v(elocity) */ 02959 /* the euler step for location then becomes */ 02960 /* x(t + dt) = x(t) + v(t~) * dt */ 02961 mul_v3_fl(dx,forcetime); 02962 02963 /* the freezer coming sooner or later */ 02964 /* 02965 if ((dot_v3v3(dx,dx)<freezeloc )&&(dot_v3v3(bp->force,bp->force)<freezeforce )){ 02966 bp->frozen /=2; 02967 } 02968 else{ 02969 bp->frozen =MIN2(bp->frozen*1.05f,1.0f); 02970 } 02971 mul_v3_fl(dx,bp->frozen); 02972 */ 02973 /* again some nasty if's to have heun in here too */ 02974 if (mode ==1){ 02975 VECCOPY(bp->prevpos,bp->pos); 02976 VECCOPY(bp->prevdx ,dx); 02977 } 02978 02979 if (mode ==2){ 02980 bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]); 02981 bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]); 02982 bp->pos[2] = bp->prevpos[2] + 0.5f * ( dx[2] + bp->prevdx[2]); 02983 maxerrpos = MAX2(maxerrpos,ABS(dx[0] - bp->prevdx[0])); 02984 maxerrpos = MAX2(maxerrpos,ABS(dx[1] - bp->prevdx[1])); 02985 maxerrpos = MAX2(maxerrpos,ABS(dx[2] - bp->prevdx[2])); 02986 02987 /* bp->choke is set when we need to pull a vertex or edge out of the collider. 02988 the collider object signals to get out by pushing hard. on the other hand 02989 we don't want to end up in deep space so we add some <viscosity> 02990 to balance that out */ 02991 if (bp->choke2 > 0.0f){ 02992 mul_v3_fl(bp->vec,(1.0f - bp->choke2)); 02993 } 02994 if (bp->choke > 0.0f){ 02995 mul_v3_fl(bp->vec,(1.0f - bp->choke)); 02996 } 02997 02998 } 02999 else { VECADD(bp->pos, bp->pos, dx);} 03000 }/*snap*/ 03001 /* so while we are looping BPs anyway do statistics on the fly */ 03002 aabbmin[0] = MIN2(aabbmin[0],bp->pos[0]); 03003 aabbmin[1] = MIN2(aabbmin[1],bp->pos[1]); 03004 aabbmin[2] = MIN2(aabbmin[2],bp->pos[2]); 03005 aabbmax[0] = MAX2(aabbmax[0],bp->pos[0]); 03006 aabbmax[1] = MAX2(aabbmax[1],bp->pos[1]); 03007 aabbmax[2] = MAX2(aabbmax[2],bp->pos[2]); 03008 if (bp->loc_flag & SBF_DOFUZZY) fuzzy =1; 03009 } /*for*/ 03010 03011 if (sb->totpoint) mul_v3_fl(cm,1.0f/sb->totpoint); 03012 if (sb->scratch){ 03013 VECCOPY(sb->scratch->aabbmin,aabbmin); 03014 VECCOPY(sb->scratch->aabbmax,aabbmax); 03015 } 03016 03017 if (err){ /* so step size will be controlled by biggest difference in slope */ 03018 if (sb->solverflags & SBSO_OLDERR) 03019 *err = MAX2(maxerrpos,maxerrvel); 03020 else 03021 *err = maxerrpos; 03022 //printf("EP %f EV %f \n",maxerrpos,maxerrvel); 03023 if (fuzzy){ 03024 *err /= sb->fuzzyness; 03025 } 03026 } 03027 } 03028 03029 /* used by heun when it overshoots */ 03030 static void softbody_restore_prev_step(Object *ob) 03031 { 03032 SoftBody *sb= ob->soft; /* is supposed to be there*/ 03033 BodyPoint *bp; 03034 int a; 03035 03036 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03037 VECCOPY(bp->vec, bp->prevvec); 03038 VECCOPY(bp->pos, bp->prevpos); 03039 } 03040 } 03041 03042 #if 0 03043 static void softbody_store_step(Object *ob) 03044 { 03045 SoftBody *sb= ob->soft; /* is supposed to be there*/ 03046 BodyPoint *bp; 03047 int a; 03048 03049 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03050 VECCOPY(bp->prevvec,bp->vec); 03051 VECCOPY(bp->prevpos,bp->pos); 03052 } 03053 } 03054 03055 03056 /* used by predictors and correctors */ 03057 static void softbody_store_state(Object *ob,float *ppos,float *pvel) 03058 { 03059 SoftBody *sb= ob->soft; /* is supposed to be there*/ 03060 BodyPoint *bp; 03061 int a; 03062 float *pp=ppos,*pv=pvel; 03063 03064 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03065 03066 VECCOPY(pv, bp->vec); 03067 pv+=3; 03068 03069 VECCOPY(pp, bp->pos); 03070 pp+=3; 03071 } 03072 } 03073 03074 /* used by predictors and correctors */ 03075 static void softbody_retrieve_state(Object *ob,float *ppos,float *pvel) 03076 { 03077 SoftBody *sb= ob->soft; /* is supposed to be there*/ 03078 BodyPoint *bp; 03079 int a; 03080 float *pp=ppos,*pv=pvel; 03081 03082 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03083 03084 VECCOPY(bp->vec,pv); 03085 pv+=3; 03086 03087 VECCOPY(bp->pos,pp); 03088 pp+=3; 03089 } 03090 } 03091 03092 /* used by predictors and correctors */ 03093 static void softbody_swap_state(Object *ob,float *ppos,float *pvel) 03094 { 03095 SoftBody *sb= ob->soft; /* is supposed to be there*/ 03096 BodyPoint *bp; 03097 int a; 03098 float *pp=ppos,*pv=pvel; 03099 float temp[3]; 03100 03101 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03102 03103 VECCOPY(temp, bp->vec); 03104 VECCOPY(bp->vec,pv); 03105 VECCOPY(pv,temp); 03106 pv+=3; 03107 03108 VECCOPY(temp, bp->pos); 03109 VECCOPY(bp->pos,pp); 03110 VECCOPY(pp,temp); 03111 pp+=3; 03112 } 03113 } 03114 #endif 03115 03116 03117 /* care for bodypoints taken out of the 'ordinary' solver step 03118 ** because they are screwed to goal by bolts 03119 ** they just need to move along with the goal in time 03120 ** we need to adjust them on sub frame timing in solver 03121 ** so now when frame is done .. put 'em to the position at the end of frame 03122 */ 03123 static void softbody_apply_goalsnap(Object *ob) 03124 { 03125 SoftBody *sb= ob->soft; /* is supposed to be there */ 03126 BodyPoint *bp; 03127 int a; 03128 03129 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03130 if (_final_goal(ob,bp) >= SOFTGOALSNAP){ 03131 VECCOPY(bp->prevpos,bp->pos); 03132 VECCOPY(bp->pos,bp->origT); 03133 } 03134 } 03135 } 03136 03137 03138 static void apply_spring_memory(Object *ob) 03139 { 03140 SoftBody *sb = ob->soft; 03141 BodySpring *bs; 03142 BodyPoint *bp1, *bp2; 03143 int a; 03144 float b,l,r; 03145 03146 if (sb && sb->totspring){ 03147 b = sb->plastic; 03148 for(a=0; a<sb->totspring; a++) { 03149 bs = &sb->bspring[a]; 03150 bp1 =&sb->bpoint[bs->v1]; 03151 bp2 =&sb->bpoint[bs->v2]; 03152 l = len_v3v3(bp1->pos,bp2->pos); 03153 r = bs->len/l; 03154 if (( r > 1.05f) || (r < 0.95)){ 03155 bs->len = ((100.0f - b) * bs->len + b*l)/100.0f; 03156 } 03157 } 03158 } 03159 } 03160 03161 /* expects full initialized softbody */ 03162 static void interpolate_exciter(Object *ob, int timescale, int time) 03163 { 03164 SoftBody *sb= ob->soft; 03165 BodyPoint *bp; 03166 float f; 03167 int a; 03168 03169 f = (float)time/(float)timescale; 03170 03171 for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { 03172 bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]); 03173 bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]); 03174 bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]); 03175 if (_final_goal(ob,bp) >= SOFTGOALSNAP){ 03176 bp->vec[0] = bp->origE[0] - bp->origS[0]; 03177 bp->vec[1] = bp->origE[1] - bp->origS[1]; 03178 bp->vec[2] = bp->origE[2] - bp->origS[2]; 03179 } 03180 } 03181 03182 } 03183 03184 03185 /* ************ convertors ********** */ 03186 03187 /* for each object type we need; 03188 - xxxx_to_softbody(Object *ob) : a full (new) copy, creates SB geometry 03189 */ 03190 03191 static void get_scalar_from_vertexgroup(Object *ob, int vertID, short groupindex, float *target) 03192 /* result 0 on success, else indicates error number 03193 -- kind of *inverse* result defintion, 03194 -- but this way we can signal error condition to caller 03195 -- and yes this function must not be here but in a *vertex group module* 03196 */ 03197 { 03198 MDeformVert *dv= NULL; 03199 int i; 03200 03201 /* spot the vert in deform vert list at mesh */ 03202 if(ob->type==OB_MESH) { 03203 Mesh *me= ob->data; 03204 if (me->dvert) 03205 dv = me->dvert + vertID; 03206 } 03207 else if(ob->type==OB_LATTICE) { /* not yet supported in softbody btw */ 03208 Lattice *lt= ob->data; 03209 if (lt->dvert) 03210 dv = lt->dvert + vertID; 03211 } 03212 if(dv) { 03213 /* Lets see if this vert is in the weight group */ 03214 for (i=0; i<dv->totweight; i++){ 03215 if (dv->dw[i].def_nr == groupindex){ 03216 *target= dv->dw[i].weight; /* got it ! */ 03217 break; 03218 } 03219 } 03220 } 03221 } 03222 03223 /* Resetting a Mesh SB object's springs */ 03224 /* Spring length are caculted from'raw' mesh vertices that are NOT altered by modifier stack. */ 03225 static void springs_from_mesh(Object *ob) 03226 { 03227 SoftBody *sb; 03228 Mesh *me= ob->data; 03229 BodyPoint *bp; 03230 int a; 03231 float scale =1.0f; 03232 03233 sb= ob->soft; 03234 if (me && sb) 03235 { 03236 /* using bp->origS as a container for spring calcualtions here 03237 ** will be overwritten sbObjectStep() to receive 03238 ** actual modifier stack positions 03239 */ 03240 if(me->totvert) { 03241 bp= ob->soft->bpoint; 03242 for(a=0; a<me->totvert; a++, bp++) { 03243 VECCOPY(bp->origS, me->mvert[a].co); 03244 mul_m4_v3(ob->obmat, bp->origS); 03245 } 03246 03247 } 03248 /* recalculate spring length for meshes here */ 03249 /* public version shrink to fit */ 03250 if (sb->springpreload != 0 ){ 03251 scale = sb->springpreload / 100.0f; 03252 } 03253 for(a=0; a<sb->totspring; a++) { 03254 BodySpring *bs = &sb->bspring[a]; 03255 bs->len= scale*len_v3v3(sb->bpoint[bs->v1].origS, sb->bpoint[bs->v2].origS); 03256 } 03257 } 03258 } 03259 03260 03261 03262 03263 /* makes totally fresh start situation */ 03264 static void mesh_to_softbody(Scene *scene, Object *ob) 03265 { 03266 SoftBody *sb; 03267 Mesh *me= ob->data; 03268 MEdge *medge= me->medge; 03269 BodyPoint *bp; 03270 BodySpring *bs; 03271 int a, totedge; 03272 if (ob->softflag & OB_SB_EDGES) totedge= me->totedge; 03273 else totedge= 0; 03274 03275 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ 03276 renew_softbody(scene, ob, me->totvert, totedge); 03277 03278 /* we always make body points */ 03279 sb= ob->soft; 03280 bp= sb->bpoint; 03281 03282 for(a=0; a<me->totvert; a++, bp++) { 03283 /* get scalar values needed *per vertex* from vertex group functions, 03284 so we can *paint* them nicly .. 03285 they are normalized [0.0..1.0] so may be we need amplitude for scale 03286 which can be done by caller but still .. i'd like it to go this way 03287 */ 03288 03289 if((ob->softflag & OB_SB_GOAL) && sb->vertgroup) { /* even this is a deprecated evil hack */ 03290 /* I'd like to have it .. if (sb->namedVG_Goal[0]) */ 03291 03292 get_scalar_from_vertexgroup(ob, a,(short) (sb->vertgroup-1), &bp->goal); 03293 /* do this always, regardless successfull read from vertex group */ 03294 /* this is where '2.5 every thing is animateable' goes wrong in the first place jow_go_for2_5 */ 03295 /* 1st coding action to take : move this to frame level */ 03296 /* reads: leave the bp->goal as it was read from vertex group / or default .. we will need it at per frame call */ 03297 /* should be fixed for meshes */ 03298 // bp->goal= sb->mingoal + bp->goal*goalfac; /* do not do here jow_go_for2_5 */ 03299 } 03300 else{ 03301 /* in consequence if no group was set .. but we want to animate it laters */ 03302 /* logically attach to goal with default first */ 03303 if(ob->softflag & OB_SB_GOAL){bp->goal =sb->defgoal;} 03304 } 03305 03306 /* to proove the concept 03307 this enables per vertex *mass painting* 03308 */ 03309 03310 if (sb->namedVG_Mass[0]) 03311 { 03312 int grp= defgroup_name_index (ob,sb->namedVG_Mass); 03313 /* printf("VGN %s %d \n",sb->namedVG_Mass,grp); */ 03314 if(grp > -1){ 03315 get_scalar_from_vertexgroup(ob, a,(short) (grp), &bp->mass); 03316 /* 2.5 bp->mass = bp->mass * sb->nodemass; */ 03317 /* printf("bp->mass %f \n",bp->mass); */ 03318 03319 } 03320 } 03321 /* first set the default */ 03322 bp->springweight = 1.0f; 03323 03324 if (sb->namedVG_Spring_K[0]) 03325 { 03326 int grp= defgroup_name_index (ob,sb->namedVG_Spring_K); 03327 //printf("VGN %s %d \n",sb->namedVG_Spring_K,grp); 03328 if(grp > -1){ 03329 get_scalar_from_vertexgroup(ob, a,(short) (grp), &bp->springweight); 03330 //printf("bp->springweight %f \n",bp->springweight); 03331 03332 } 03333 } 03334 03335 03336 } 03337 03338 /* but we only optionally add body edge springs */ 03339 if (ob->softflag & OB_SB_EDGES) { 03340 if(medge) { 03341 bs= sb->bspring; 03342 for(a=me->totedge; a>0; a--, medge++, bs++) { 03343 bs->v1= medge->v1; 03344 bs->v2= medge->v2; 03345 bs->springtype=SB_EDGE; 03346 } 03347 03348 03349 /* insert *diagonal* springs in quads if desired */ 03350 if (ob->softflag & OB_SB_QUADS) { 03351 add_mesh_quad_diag_springs(ob); 03352 } 03353 03354 build_bps_springlist(ob); /* scan for springs attached to bodypoints ONCE */ 03355 /* insert *other second order* springs if desired */ 03356 if (sb->secondspring > 0.0000001f) { 03357 add_2nd_order_springs(ob,sb->secondspring); /* exploits the the first run of build_bps_springlist(ob);*/ 03358 build_bps_springlist(ob); /* yes we need to do it again*/ 03359 } 03360 springs_from_mesh(ob); /* write the 'rest'-length of the springs */ 03361 if (ob->softflag & OB_SB_SELF) {calculate_collision_balls(ob);} 03362 03363 } 03364 03365 } 03366 03367 } 03368 03369 static void mesh_faces_to_scratch(Object *ob) 03370 { 03371 SoftBody *sb= ob->soft; 03372 Mesh *me= ob->data; 03373 MFace *mface; 03374 BodyFace *bodyface; 03375 int a; 03376 /* alloc and copy faces*/ 03377 03378 bodyface = sb->scratch->bodyface = MEM_mallocN(sizeof(BodyFace)*me->totface,"SB_body_Faces"); 03379 //memcpy(sb->scratch->mface,me->mface,sizeof(MFace)*me->totface); 03380 mface = me->mface; 03381 for(a=0; a<me->totface; a++, mface++, bodyface++) { 03382 bodyface->v1 = mface->v1; 03383 bodyface->v2 = mface->v2; 03384 bodyface->v3 = mface->v3; 03385 bodyface->v4 = mface->v4; 03386 bodyface->ext_force[0] = bodyface->ext_force[1] = bodyface->ext_force[2] = 0.0f; 03387 bodyface->flag =0; 03388 } 03389 sb->scratch->totface = me->totface; 03390 } 03391 static void reference_to_scratch(Object *ob) 03392 { 03393 SoftBody *sb= ob->soft; 03394 ReferenceVert *rp; 03395 BodyPoint *bp; 03396 float accu_pos[3] ={0.f,0.f,0.f}; 03397 float accu_mass = 0.f; 03398 int a; 03399 03400 sb->scratch->Ref.ivert = MEM_mallocN(sizeof(ReferenceVert)*sb->totpoint,"SB_Reference"); 03401 bp= ob->soft->bpoint; 03402 rp= sb->scratch->Ref.ivert; 03403 for(a=0; a<sb->totpoint; a++, rp++, bp++) { 03404 VECCOPY(rp->pos,bp->pos); 03405 VECADD(accu_pos,accu_pos,bp->pos); 03406 accu_mass += _final_mass(ob,bp); 03407 } 03408 mul_v3_fl(accu_pos,1.0f/accu_mass); 03409 VECCOPY(sb->scratch->Ref.com,accu_pos); 03410 /* printf("reference_to_scratch \n"); */ 03411 } 03412 03413 /* 03414 helper function to get proper spring length 03415 when object is rescaled 03416 */ 03417 static float globallen(float *v1,float *v2,Object *ob) 03418 { 03419 float p1[3],p2[3]; 03420 VECCOPY(p1,v1); 03421 mul_m4_v3(ob->obmat, p1); 03422 VECCOPY(p2,v2); 03423 mul_m4_v3(ob->obmat, p2); 03424 return len_v3v3(p1,p2); 03425 } 03426 03427 static void makelatticesprings(Lattice *lt, BodySpring *bs, int dostiff,Object *ob) 03428 { 03429 BPoint *bp=lt->def, *bpu; 03430 int u, v, w, dv, dw, bpc=0, bpuc; 03431 03432 dv= lt->pntsu; 03433 dw= dv*lt->pntsv; 03434 03435 for(w=0; w<lt->pntsw; w++) { 03436 03437 for(v=0; v<lt->pntsv; v++) { 03438 03439 for(u=0, bpuc=0, bpu=NULL; u<lt->pntsu; u++, bp++, bpc++) { 03440 03441 if(w) { 03442 bs->v1 = bpc; 03443 bs->v2 = bpc-dw; 03444 bs->springtype=SB_EDGE; 03445 bs->len= globallen((bp-dw)->vec, bp->vec,ob); 03446 bs++; 03447 } 03448 if(v) { 03449 bs->v1 = bpc; 03450 bs->v2 = bpc-dv; 03451 bs->springtype=SB_EDGE; 03452 bs->len= globallen((bp-dv)->vec, bp->vec,ob); 03453 bs++; 03454 } 03455 if(u) { 03456 bs->v1 = bpuc; 03457 bs->v2 = bpc; 03458 bs->springtype=SB_EDGE; 03459 bs->len= globallen((bpu)->vec, bp->vec,ob); 03460 bs++; 03461 } 03462 03463 if (dostiff) { 03464 03465 if(w){ 03466 if( v && u ) { 03467 bs->v1 = bpc; 03468 bs->v2 = bpc-dw-dv-1; 03469 bs->springtype=SB_BEND; 03470 bs->len= globallen((bp-dw-dv-1)->vec, bp->vec,ob); 03471 bs++; 03472 } 03473 if( (v < lt->pntsv-1) && (u) ) { 03474 bs->v1 = bpc; 03475 bs->v2 = bpc-dw+dv-1; 03476 bs->springtype=SB_BEND; 03477 bs->len= globallen((bp-dw+dv-1)->vec, bp->vec,ob); 03478 bs++; 03479 } 03480 } 03481 03482 if(w < lt->pntsw -1){ 03483 if( v && u ) { 03484 bs->v1 = bpc; 03485 bs->v2 = bpc+dw-dv-1; 03486 bs->springtype=SB_BEND; 03487 bs->len= globallen((bp+dw-dv-1)->vec, bp->vec,ob); 03488 bs++; 03489 } 03490 if( (v < lt->pntsv-1) && (u) ) { 03491 bs->v1 = bpc; 03492 bs->v2 = bpc+dw+dv-1; 03493 bs->springtype=SB_BEND; 03494 bs->len= globallen((bp+dw+dv-1)->vec, bp->vec,ob); 03495 bs++; 03496 } 03497 } 03498 } 03499 bpu = bp; 03500 bpuc = bpc; 03501 } 03502 } 03503 } 03504 } 03505 03506 03507 /* makes totally fresh start situation */ 03508 static void lattice_to_softbody(Scene *scene, Object *ob) 03509 { 03510 Lattice *lt= ob->data; 03511 SoftBody *sb; 03512 int totvert, totspring = 0; 03513 03514 totvert= lt->pntsu*lt->pntsv*lt->pntsw; 03515 03516 if (ob->softflag & OB_SB_EDGES){ 03517 totspring = ((lt->pntsu -1) * lt->pntsv 03518 + (lt->pntsv -1) * lt->pntsu) * lt->pntsw 03519 +lt->pntsu*lt->pntsv*(lt->pntsw -1); 03520 if (ob->softflag & OB_SB_QUADS){ 03521 totspring += 4*(lt->pntsu -1) * (lt->pntsv -1) * (lt->pntsw-1); 03522 } 03523 } 03524 03525 03526 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ 03527 renew_softbody(scene, ob, totvert, totspring); 03528 sb= ob->soft; /* can be created in renew_softbody() */ 03529 03530 /* weights from bpoints, same code used as for mesh vertices */ 03531 /* if((ob->softflag & OB_SB_GOAL) && sb->vertgroup) { 2.4x one*/ 03532 /* new! take the weights from lattice vertex anyhow */ 03533 if(ob->softflag & OB_SB_GOAL){ 03534 BodyPoint *bp= sb->bpoint; 03535 BPoint *bpnt= lt->def; 03536 /* jow_go_for2_5 */ 03537 int a; 03538 03539 for(a=0; a<totvert; a++, bp++, bpnt++) { 03540 bp->goal= bpnt->weight; 03541 } 03542 } 03543 03544 /* create some helper edges to enable SB lattice to be usefull at all */ 03545 if (ob->softflag & OB_SB_EDGES){ 03546 makelatticesprings(lt,ob->soft->bspring,ob->softflag & OB_SB_QUADS,ob); 03547 build_bps_springlist(ob); /* link bps to springs */ 03548 } 03549 } 03550 03551 /* makes totally fresh start situation */ 03552 static void curve_surf_to_softbody(Scene *scene, Object *ob) 03553 { 03554 Curve *cu= ob->data; 03555 SoftBody *sb; 03556 BodyPoint *bp; 03557 BodySpring *bs; 03558 Nurb *nu; 03559 BezTriple *bezt; 03560 BPoint *bpnt; 03561 int a, curindex=0; 03562 int totvert, totspring = 0, setgoal=0; 03563 03564 totvert= count_curveverts(&cu->nurb); 03565 03566 if (ob->softflag & OB_SB_EDGES){ 03567 if(ob->type==OB_CURVE) { 03568 totspring= totvert - BLI_countlist(&cu->nurb); 03569 } 03570 } 03571 03572 /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ 03573 renew_softbody(scene, ob, totvert, totspring); 03574 sb= ob->soft; /* can be created in renew_softbody() */ 03575 03576 /* set vars now */ 03577 bp= sb->bpoint; 03578 bs= sb->bspring; 03579 03580 /* weights from bpoints, same code used as for mesh vertices */ 03581 /* if((ob->softflag & OB_SB_GOAL) && sb->vertgroup) 2.4x hack*/ 03582 /* new! take the weights from curve vertex anyhow */ 03583 if(ob->softflag & OB_SB_GOAL) 03584 setgoal= 1; 03585 03586 for(nu= cu->nurb.first; nu; nu= nu->next) { 03587 if(nu->bezt) { 03588 /* bezier case ; this is nicly said naive; who ever wrote this part, it was not me (JOW) :) */ 03589 /* a: never ever make tangent handles (sub) and or (ob)ject to collision */ 03590 /* b: rather calculate them using some C2 (C2= continous in second derivate -> no jump in bending ) condition */ 03591 /* not too hard to do, but needs some more code to care for; some one may want look at it JOW 2010/06/12*/ 03592 for(bezt=nu->bezt, a=0; a<nu->pntsu; a++, bezt++, bp+=3, curindex+=3) { 03593 if(setgoal) { 03594 bp->goal= bezt->weight; 03595 03596 /* all three triples */ 03597 (bp+1)->goal= bp->goal; 03598 (bp+2)->goal= bp->goal; 03599 /*do not collide handles */ 03600 (bp+1)->loc_flag |= SBF_OUTOFCOLLISION; 03601 (bp+2)->loc_flag |= SBF_OUTOFCOLLISION; 03602 } 03603 03604 if(totspring) { 03605 if(a>0) { 03606 bs->v1= curindex-3; 03607 bs->v2= curindex; 03608 bs->springtype=SB_HANDLE; 03609 bs->len= globallen( (bezt-1)->vec[0], bezt->vec[0], ob ); 03610 bs++; 03611 } 03612 bs->v1= curindex; 03613 bs->v2= curindex+1; 03614 bs->springtype=SB_HANDLE; 03615 bs->len= globallen( bezt->vec[0], bezt->vec[1], ob ); 03616 bs++; 03617 03618 bs->v1= curindex+1; 03619 bs->v2= curindex+2; 03620 bs->springtype=SB_HANDLE; 03621 bs->len= globallen( bezt->vec[1], bezt->vec[2], ob ); 03622 bs++; 03623 } 03624 } 03625 } 03626 else { 03627 for(bpnt=nu->bp, a=0; a<nu->pntsu*nu->pntsv; a++, bpnt++, bp++, curindex++) { 03628 if(setgoal) { 03629 bp->goal= bpnt->weight; 03630 } 03631 if(totspring && a>0) { 03632 bs->v1= curindex-1; 03633 bs->v2= curindex; 03634 bs->springtype=SB_EDGE; 03635 bs->len= globallen( (bpnt-1)->vec, bpnt->vec , ob ); 03636 bs++; 03637 } 03638 } 03639 } 03640 } 03641 03642 if(totspring) 03643 { 03644 build_bps_springlist(ob); /* link bps to springs */ 03645 if (ob->softflag & OB_SB_SELF) {calculate_collision_balls(ob);} 03646 } 03647 } 03648 03649 /* copies softbody result back in object */ 03650 static void softbody_to_object(Object *ob, float (*vertexCos)[3], int numVerts, int local) 03651 { 03652 SoftBody *sb= ob->soft; 03653 if(sb){ 03654 BodyPoint *bp= sb->bpoint; 03655 int a; 03656 if(sb->solverflags & SBSO_ESTIMATEIPO){SB_estimate_transform(ob,sb->lcom,sb->lrot,sb->lscale);} 03657 /* inverse matrix is not uptodate... */ 03658 invert_m4_m4(ob->imat, ob->obmat); 03659 03660 for(a=0; a<numVerts; a++, bp++) { 03661 VECCOPY(vertexCos[a], bp->pos); 03662 if(local==0) 03663 mul_m4_v3(ob->imat, vertexCos[a]); /* softbody is in global coords, baked optionally not */ 03664 } 03665 } 03666 } 03667 03668 /* +++ ************ maintaining scratch *************** */ 03669 static void sb_new_scratch(SoftBody *sb) 03670 { 03671 if (!sb) return; 03672 sb->scratch = MEM_callocN(sizeof(SBScratch), "SBScratch"); 03673 sb->scratch->colliderhash = BLI_ghash_new(BLI_ghashutil_ptrhash,BLI_ghashutil_ptrcmp, "sb_new_scratch gh"); 03674 sb->scratch->bodyface = NULL; 03675 sb->scratch->totface = 0; 03676 sb->scratch->aabbmax[0]=sb->scratch->aabbmax[1]=sb->scratch->aabbmax[2] = 1.0e30f; 03677 sb->scratch->aabbmin[0]=sb->scratch->aabbmin[1]=sb->scratch->aabbmin[2] = -1.0e30f; 03678 sb->scratch->Ref.ivert = NULL; 03679 03680 } 03681 /* --- ************ maintaining scratch *************** */ 03682 03683 /* ************ Object level, exported functions *************** */ 03684 03685 /* allocates and initializes general main data */ 03686 SoftBody *sbNew(Scene *scene) 03687 { 03688 SoftBody *sb; 03689 03690 sb= MEM_callocN(sizeof(SoftBody), "softbody"); 03691 03692 sb->mediafrict= 0.5f; 03693 sb->nodemass= 1.0f; 03694 sb->grav= 9.8f; 03695 sb->physics_speed= 1.0f; 03696 sb->rklimit= 0.1f; 03697 03698 sb->goalspring= 0.5f; 03699 sb->goalfrict= 0.0f; 03700 sb->mingoal= 0.0f; 03701 sb->maxgoal= 1.0f; 03702 sb->defgoal= 0.7f; 03703 03704 sb->inspring= 0.5f; 03705 sb->infrict= 0.5f; 03706 /*todo backward file compat should copy inspring to inpush while reading old files*/ 03707 sb->inpush = 0.5f; 03708 03709 sb->interval= 10; 03710 sb->sfra= scene->r.sfra; 03711 sb->efra= scene->r.efra; 03712 03713 sb->colball = 0.49f; 03714 sb->balldamp = 0.50f; 03715 sb->ballstiff= 1.0f; 03716 sb->sbc_mode = 1; 03717 03718 03719 sb->minloops = 10; 03720 sb->maxloops = 300; 03721 03722 sb->choke = 3; 03723 sb_new_scratch(sb); 03724 /*todo backward file compat should set sb->shearstiff = 1.0f while reading old files*/ 03725 sb->shearstiff = 1.0f; 03726 sb->solverflags |= SBSO_OLDERR; 03727 03728 sb->pointcache = BKE_ptcache_add(&sb->ptcaches); 03729 03730 if(!sb->effector_weights) 03731 sb->effector_weights = BKE_add_effector_weights(NULL); 03732 03733 return sb; 03734 } 03735 03736 /* frees all */ 03737 void sbFree(SoftBody *sb) 03738 { 03739 free_softbody_intern(sb); 03740 BKE_ptcache_free_list(&sb->ptcaches); 03741 sb->pointcache = NULL; 03742 if(sb->effector_weights) 03743 MEM_freeN(sb->effector_weights); 03744 MEM_freeN(sb); 03745 } 03746 03747 void sbFreeSimulation(SoftBody *sb) 03748 { 03749 free_softbody_intern(sb); 03750 } 03751 03752 /* makes totally fresh start situation */ 03753 void sbObjectToSoftbody(Object *ob) 03754 { 03755 //ob->softflag |= OB_SB_REDO; 03756 03757 free_softbody_intern(ob->soft); 03758 } 03759 03760 static int object_has_edges(Object *ob) 03761 { 03762 if(ob->type==OB_MESH) { 03763 return ((Mesh*) ob->data)->totedge; 03764 } 03765 else if(ob->type==OB_LATTICE) { 03766 return 1; 03767 } 03768 else { 03769 return 0; 03770 } 03771 } 03772 03773 /* SB global visible functions */ 03774 void sbSetInterruptCallBack(int (*f)(void)) 03775 { 03776 SB_localInterruptCallBack = f; 03777 } 03778 03779 static void softbody_update_positions(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) 03780 { 03781 BodyPoint *bp; 03782 int a; 03783 03784 if(!sb || !sb->bpoint) 03785 return; 03786 03787 for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { 03788 /* store where goals are now */ 03789 VECCOPY(bp->origS, bp->origE); 03790 /* copy the position of the goals at desired end time */ 03791 VECCOPY(bp->origE, vertexCos[a]); 03792 /* vertexCos came from local world, go global */ 03793 mul_m4_v3(ob->obmat, bp->origE); 03794 /* just to be save give bp->origT a defined value 03795 will be calulated in interpolate_exciter()*/ 03796 VECCOPY(bp->origT, bp->origE); 03797 } 03798 } 03799 03800 03801 /* void SB_estimate_transform */ 03802 /* input Object *ob out (says any object that can do SB like mesh,lattice,curve ) 03803 output float lloc[3],float lrot[3][3],float lscale[3][3] 03804 that is: 03805 a precise position vector denoting the motion of the center of mass 03806 give a rotation/scale matrix using averaging method, that's why estimate and not calculate 03807 see: this is kind of reverse engeneering: having to states of a point cloud and recover what happend 03808 our advantage here we know the identity of the vertex 03809 there are others methods giving other results. 03810 lloc,lrot,lscale are allowed to be NULL, just in case you don't need it. 03811 should be pretty useful for pythoneers :) 03812 not! velocity .. 2nd order stuff 03813 vcloud_estimate_transform see 03814 */ 03815 03816 void SB_estimate_transform(Object *ob,float lloc[3],float lrot[3][3],float lscale[3][3]) 03817 { 03818 BodyPoint *bp; 03819 ReferenceVert *rp; 03820 SoftBody *sb = NULL; 03821 float (*opos)[3]; 03822 float (*rpos)[3]; 03823 float com[3],rcom[3]; 03824 int a; 03825 03826 if(!ob ||!ob->soft) return; /* why did we get here ? */ 03827 sb= ob->soft; 03828 if(!sb || !sb->bpoint) return; 03829 opos= MEM_callocN( (sb->totpoint)*3*sizeof(float), "SB_OPOS"); 03830 rpos= MEM_callocN( (sb->totpoint)*3*sizeof(float), "SB_RPOS"); 03831 /* might filter vertex selection with a vertex group */ 03832 for(a=0, bp=sb->bpoint, rp=sb->scratch->Ref.ivert; a<sb->totpoint; a++, bp++, rp++) { 03833 VECCOPY(rpos[a],rp->pos); 03834 VECCOPY(opos[a],bp->pos); 03835 } 03836 03837 vcloud_estimate_transform(sb->totpoint, opos, NULL, rpos, NULL, com, rcom,lrot,lscale); 03838 //VECSUB(com,com,rcom); 03839 if (lloc) VECCOPY(lloc,com); 03840 VECCOPY(sb->lcom,com); 03841 if (lscale) copy_m3_m3(sb->lscale,lscale); 03842 if (lrot) copy_m3_m3(sb->lrot,lrot); 03843 03844 03845 MEM_freeN(opos); 03846 MEM_freeN(rpos); 03847 } 03848 03849 static void softbody_reset(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) 03850 { 03851 BodyPoint *bp; 03852 int a; 03853 03854 for(a=0,bp=sb->bpoint; a<numVerts; a++, bp++) { 03855 VECCOPY(bp->pos, vertexCos[a]); 03856 mul_m4_v3(ob->obmat, bp->pos); /* yep, sofbody is global coords*/ 03857 VECCOPY(bp->origS, bp->pos); 03858 VECCOPY(bp->origE, bp->pos); 03859 VECCOPY(bp->origT, bp->pos); 03860 bp->vec[0]= bp->vec[1]= bp->vec[2]= 0.0f; 03861 03862 /* the bp->prev*'s are for rolling back from a canceled try to propagate in time 03863 adaptive step size algo in a nutshell: 03864 1. set sheduled time step to new dtime 03865 2. try to advance the sheduled time step, beeing optimistic execute it 03866 3. check for success 03867 3.a we 're fine continue, may be we can increase sheduled time again ?? if so, do so! 03868 3.b we did exceed error limit --> roll back, shorten the sheduled time and try again at 2. 03869 4. check if we did reach dtime 03870 4.a nope we need to do some more at 2. 03871 4.b yup we're done 03872 */ 03873 03874 VECCOPY(bp->prevpos, bp->pos); 03875 VECCOPY(bp->prevvec, bp->vec); 03876 VECCOPY(bp->prevdx, bp->vec); 03877 VECCOPY(bp->prevdv, bp->vec); 03878 } 03879 03880 /* make a nice clean scratch struc */ 03881 free_scratch(sb); /* clear if any */ 03882 sb_new_scratch(sb); /* make a new */ 03883 sb->scratch->needstobuildcollider=1; 03884 zero_v3(sb->lcom); 03885 unit_m3(sb->lrot); 03886 unit_m3(sb->lscale); 03887 03888 03889 03890 /* copy some info to scratch */ 03891 /* we only need that if we want to reconstruct IPO */ 03892 if(1) { 03893 reference_to_scratch(ob); 03894 SB_estimate_transform(ob,NULL,NULL,NULL); 03895 SB_estimate_transform(ob,NULL,NULL,NULL); 03896 } 03897 switch(ob->type) { 03898 case OB_MESH: 03899 if (ob->softflag & OB_SB_FACECOLL) mesh_faces_to_scratch(ob); 03900 break; 03901 case OB_LATTICE: 03902 break; 03903 case OB_CURVE: 03904 case OB_SURF: 03905 break; 03906 default: 03907 break; 03908 } 03909 } 03910 03911 static void softbody_step(Scene *scene, Object *ob, SoftBody *sb, float dtime) 03912 { 03913 /* the simulator */ 03914 float forcetime; 03915 double sct,sst; 03916 03917 03918 sst=PIL_check_seconds_timer(); 03919 /* Integration back in time is possible in theory, but pretty useless here. 03920 So we refuse to do so. Since we do not know anything about 'outside' canges 03921 especially colliders we refuse to go more than 10 frames. 03922 */ 03923 if(dtime < 0 || dtime > 10.5f) return; 03924 03925 ccd_update_deflector_hash(scene, ob, sb->scratch->colliderhash); 03926 03927 if(sb->scratch->needstobuildcollider){ 03928 if (query_external_colliders(scene, ob)){ 03929 ccd_build_deflector_hash(scene, ob, sb->scratch->colliderhash); 03930 } 03931 sb->scratch->needstobuildcollider=0; 03932 } 03933 03934 if (sb->solver_ID < 2) { 03935 /* special case of 2nd order Runge-Kutta type AKA Heun */ 03936 int mid_flags=0; 03937 float err = 0; 03938 float forcetimemax = 1.0f; /* set defaults guess we shall do one frame */ 03939 float forcetimemin = 0.01f; /* set defaults guess 1/100 is tight enough */ 03940 float timedone =0.0; /* how far did we get without violating error condition */ 03941 /* loops = counter for emergency brake 03942 * we don't want to lock up the system if physics fail 03943 */ 03944 int loops =0 ; 03945 03946 SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */ 03947 /* adjust loop limits */ 03948 if (sb->minloops > 0) forcetimemax = dtime / sb->minloops; 03949 if (sb->maxloops > 0) forcetimemin = dtime / sb->maxloops; 03950 03951 if(sb->solver_ID>0) mid_flags |= MID_PRESERVE; 03952 03953 forcetime =forcetimemax; /* hope for integrating in one step */ 03954 while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ) 03955 { 03956 /* set goals in time */ 03957 interpolate_exciter(ob,200,(int)(200.0*(timedone/dtime))); 03958 03959 sb->scratch->flag &= ~SBF_DOFUZZY; 03960 /* do predictive euler step */ 03961 softbody_calc_forces(scene, ob, forcetime,timedone/dtime,0); 03962 03963 softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags); 03964 03965 /* crop new slope values to do averaged slope step */ 03966 softbody_calc_forces(scene, ob, forcetime,timedone/dtime,0); 03967 03968 softbody_apply_forces(ob, forcetime, 2, &err,mid_flags); 03969 softbody_apply_goalsnap(ob); 03970 03971 if (err > SoftHeunTol) { /* error needs to be scaled to some quantity */ 03972 03973 if (forcetime > forcetimemin){ 03974 forcetime = MAX2(forcetime / 2.0f,forcetimemin); 03975 softbody_restore_prev_step(ob); 03976 //printf("down,"); 03977 } 03978 else { 03979 timedone += forcetime; 03980 } 03981 } 03982 else { 03983 float newtime = forcetime * 1.1f; /* hope for 1.1 times better conditions in next step */ 03984 03985 if (sb->scratch->flag & SBF_DOFUZZY){ 03986 //if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */ 03987 newtime = forcetime; 03988 //} 03989 } 03990 else { 03991 if (err > SoftHeunTol/2.0f) { /* stay with this stepsize unless err really small */ 03992 newtime = forcetime; 03993 } 03994 } 03995 timedone += forcetime; 03996 newtime=MIN2(forcetimemax,MAX2(newtime,forcetimemin)); 03997 //if (newtime > forcetime) printf("up,"); 03998 if (forcetime > 0.0) 03999 forcetime = MIN2(dtime - timedone,newtime); 04000 else 04001 forcetime = MAX2(dtime - timedone,newtime); 04002 } 04003 loops++; 04004 if(sb->solverflags & SBSO_MONITOR ){ 04005 sct=PIL_check_seconds_timer(); 04006 if (sct-sst > 0.5f) printf("%3.0f%% \r",100.0f*timedone/dtime); 04007 } 04008 /* ask for user break */ 04009 if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break; 04010 04011 } 04012 /* move snapped to final position */ 04013 interpolate_exciter(ob, 2, 2); 04014 softbody_apply_goalsnap(ob); 04015 04016 // if(G.f & G_DEBUG){ 04017 if(sb->solverflags & SBSO_MONITOR ){ 04018 if (loops > HEUNWARNLIMIT) /* monitor high loop counts */ 04019 printf("\r needed %d steps/frame",loops); 04020 } 04021 04022 } 04023 else if (sb->solver_ID == 2) 04024 {/* do semi "fake" implicit euler */ 04025 //removed 04026 }/*SOLVER SELECT*/ 04027 else if (sb->solver_ID == 4) 04028 { 04029 /* do semi "fake" implicit euler */ 04030 }/*SOLVER SELECT*/ 04031 else if (sb->solver_ID == 3){ 04032 /* do "stupid" semi "fake" implicit euler */ 04033 //removed 04034 04035 }/*SOLVER SELECT*/ 04036 else{ 04037 printf("softbody no valid solver ID!"); 04038 }/*SOLVER SELECT*/ 04039 if(sb->plastic){ apply_spring_memory(ob);} 04040 04041 if(sb->solverflags & SBSO_MONITOR ){ 04042 sct=PIL_check_seconds_timer(); 04043 if ((sct-sst > 0.5f) || (G.f & G_DEBUG)) printf(" solver time %f sec %s \n",sct-sst,ob->id.name); 04044 } 04045 } 04046 04047 /* simulates one step. framenr is in frames */ 04048 void sbObjectStep(Scene *scene, Object *ob, float cfra, float (*vertexCos)[3], int numVerts) 04049 { 04050 SoftBody *sb= ob->soft; 04051 PointCache *cache; 04052 PTCacheID pid; 04053 float dtime, timescale; 04054 int framedelta, framenr, startframe, endframe; 04055 int cache_result; 04056 04057 cache= sb->pointcache; 04058 04059 framenr= (int)cfra; 04060 framedelta= framenr - cache->simframe; 04061 04062 BKE_ptcache_id_from_softbody(&pid, ob, sb); 04063 BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, ×cale); 04064 04065 /* check for changes in mesh, should only happen in case the mesh 04066 * structure changes during an animation */ 04067 if(sb->bpoint && numVerts != sb->totpoint) { 04068 BKE_ptcache_invalidate(cache); 04069 return; 04070 } 04071 04072 /* clamp frame ranges */ 04073 if(framenr < startframe) { 04074 BKE_ptcache_invalidate(cache); 04075 return; 04076 } 04077 else if(framenr > endframe) { 04078 framenr = endframe; 04079 } 04080 04081 /* verify if we need to create the softbody data */ 04082 if(sb->bpoint == NULL || 04083 ((ob->softflag & OB_SB_EDGES) && !ob->soft->bspring && object_has_edges(ob))) { 04084 04085 switch(ob->type) { 04086 case OB_MESH: 04087 mesh_to_softbody(scene, ob); 04088 break; 04089 case OB_LATTICE: 04090 lattice_to_softbody(scene, ob); 04091 break; 04092 case OB_CURVE: 04093 case OB_SURF: 04094 curve_surf_to_softbody(scene, ob); 04095 break; 04096 default: 04097 renew_softbody(scene, ob, numVerts, 0); 04098 break; 04099 } 04100 04101 softbody_update_positions(ob, sb, vertexCos, numVerts); 04102 softbody_reset(ob, sb, vertexCos, numVerts); 04103 } 04104 04105 /* continue physics special case */ 04106 if(BKE_ptcache_get_continue_physics()) { 04107 BKE_ptcache_invalidate(cache); 04108 /* do simulation */ 04109 dtime = timescale; 04110 softbody_update_positions(ob, sb, vertexCos, numVerts); 04111 softbody_step(scene, ob, sb, dtime); 04112 softbody_to_object(ob, vertexCos, numVerts, 0); 04113 return; 04114 } 04115 04116 /* still no points? go away */ 04117 if(sb->totpoint==0) { 04118 return; 04119 } 04120 if(framenr == startframe) { 04121 BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED); 04122 04123 /* first frame, no simulation to do, just set the positions */ 04124 softbody_update_positions(ob, sb, vertexCos, numVerts); 04125 04126 BKE_ptcache_validate(cache, framenr); 04127 cache->flag &= ~PTCACHE_REDO_NEEDED; 04128 return; 04129 } 04130 04131 /* try to read from cache */ 04132 cache_result = BKE_ptcache_read(&pid, framenr); 04133 04134 if(cache_result == PTCACHE_READ_EXACT || cache_result == PTCACHE_READ_INTERPOLATED) { 04135 softbody_to_object(ob, vertexCos, numVerts, sb->local); 04136 04137 BKE_ptcache_validate(cache, framenr); 04138 04139 if(cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED) 04140 BKE_ptcache_write(&pid, framenr); 04141 04142 return; 04143 } 04144 else if(cache_result==PTCACHE_READ_OLD) { 04145 ; /* do nothing */ 04146 } 04147 else if(/*ob->id.lib || */(cache->flag & PTCACHE_BAKED)) { /* "library linking & pointcaches" has to be solved properly at some point */ 04148 /* if baked and nothing in cache, do nothing */ 04149 BKE_ptcache_invalidate(cache); 04150 return; 04151 } 04152 04153 /* if on second frame, write cache for first frame */ 04154 if(cache->simframe == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact==0)) 04155 BKE_ptcache_write(&pid, startframe); 04156 04157 softbody_update_positions(ob, sb, vertexCos, numVerts); 04158 04159 /* checking time: */ 04160 dtime = framedelta*timescale; 04161 04162 /* do simulation */ 04163 softbody_step(scene, ob, sb, dtime); 04164 04165 softbody_to_object(ob, vertexCos, numVerts, 0); 04166 04167 BKE_ptcache_validate(cache, framenr); 04168 BKE_ptcache_write(&pid, framenr); 04169 } 04170