/* * ***** BEGIN GPL LICENSE BLOCK ***** * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software Foundation, * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * * The Original Code is Copyright (C) Blender Foundation * All rights reserved. * * The Original Code is: all of this file. * * Contributor(s): none yet. * * ***** END GPL LICENSE BLOCK ***** */ /** \file blender/blenkernel/intern/softbody.c * \ingroup bke */ /* ****** variables on the UI for now float mediafrict; friction to env float nodemass; softbody mass of *vertex* float grav; softbody amount of gravitaion to apply float goalspring; softbody goal springs float goalfrict; softbody goal springs friction float mingoal; quick limits for goal float maxgoal; float inspring; softbody inner springs float infrict; softbody inner springs friction ***** */ #include #include #include #include "MEM_guardedalloc.h" /* types */ #include "DNA_object_types.h" #include "DNA_scene_types.h" #include "DNA_lattice_types.h" #include "DNA_curve_types.h" #include "DNA_mesh_types.h" #include "DNA_meshdata_types.h" #include "BLI_math.h" #include "BLI_utildefines.h" #include "BLI_listbase.h" #include "BLI_ghash.h" #include "BLI_threads.h" #include "BKE_curve.h" #include "BKE_effect.h" #include "BKE_global.h" #include "BKE_modifier.h" #include "BKE_softbody.h" #include "BKE_pointcache.h" #include "BKE_deform.h" #include "BKE_mesh.h" #include "BKE_scene.h" #include "PIL_time.h" // #include "ONL_opennl.h" remove linking to ONL for now /* callbacks for errors and interrupts and some goo */ static int (*SB_localInterruptCallBack)(void) = NULL; /* ********** soft body engine ******* */ typedef enum {SB_EDGE=1, SB_BEND=2, SB_STIFFQUAD=3, SB_HANDLE=4} type_spring; typedef struct BodySpring { int v1, v2; float len, cf, load; float ext_force[3]; /* edges colliding and sailing */ type_spring springtype; short flag; } BodySpring; typedef struct BodyFace { int v1, v2, v3, v4; float ext_force[3]; /* faces colliding */ short flag; } BodyFace; typedef struct ReferenceVert { float pos[3]; /* position relative to com */ float mass; /* node mass */ } ReferenceVert; typedef struct ReferenceState { float com[3]; /* center of mass*/ ReferenceVert *ivert; /* list of intial values */ } ReferenceState; /*private scratch pad for caching and other data only needed when alive*/ typedef struct SBScratch { GHash *colliderhash; short needstobuildcollider; short flag; BodyFace *bodyface; int totface; float aabbmin[3], aabbmax[3]; ReferenceState Ref; } SBScratch; typedef struct SB_thread_context { Scene *scene; Object *ob; float forcetime; float timenow; int ifirst; int ilast; ListBase *do_effector; int do_deflector; float fieldfactor; float windfactor; int nr; int tot; } SB_thread_context; #define NLF_BUILD 1 #if 0 # define NLF_SOLVE 2 #endif #define MID_PRESERVE 1 #define SOFTGOALSNAP 0.999f /* if bp-> goal is above make it a *forced follow original* and skip all ODE stuff for this bp * removes *unnecessary* stiffness from ODE system */ #define HEUNWARNLIMIT 1 /* 500 would be fine i think for detecting severe *stiff* stuff */ #define BSF_INTERSECT 1 /* edge intersects collider face */ /* private definitions for bodypoint states */ #define SBF_DOFUZZY 1 /* Bodypoint do fuzzy */ #define SBF_OUTOFCOLLISION 2 /* Bodypoint does not collide */ #define BFF_INTERSECT 1 /* collider edge intrudes face */ #define BFF_CLOSEVERT 2 /* collider vertex repulses face */ static float SoftHeunTol = 1.0f; /* humm .. this should be calculated from sb parameters and sizes */ /* local prototypes */ static void free_softbody_intern(SoftBody *sb); /* aye this belongs to arith.c */ static void Vec3PlusStVec(float *v, float s, float *v1); /*+++ frame based timing +++*/ /*physical unit of force is [kg * m / sec^2]*/ static float sb_grav_force_scale(Object *UNUSED(ob)) /* since unit of g is [m/sec^2] and F = mass * g we rescale unit mass of node to 1 gramm * put it to a function here, so we can add user options later without touching simulation code */ { return (0.001f); } static float sb_fric_force_scale(Object *UNUSED(ob)) /* rescaling unit of drag [1 / sec] to somehow reasonable * put it to a function here, so we can add user options later without touching simulation code */ { return (0.01f); } static float sb_time_scale(Object *ob) /* defining the frames to *real* time relation */ { SoftBody *sb= ob->soft; /* is supposed to be there */ if (sb) { return(sb->physics_speed); /* hrms .. this could be IPO as well :) * estimated range [0.001 sluggish slug - 100.0 very fast (i hope ODE solver can handle that)] * 1 approx = a unit 1 pendulum at g = 9.8 [earth conditions] has period 65 frames * theory would give a 50 frames period .. so there must be something inaccurate .. looking for that (BM) */ } return (1.0f); /* * this would be frames/sec independent timing assuming 25 fps is default * but does not work very well with NLA * return (25.0f/scene->r.frs_sec) */ } /*--- frame based timing ---*/ /* helper functions for everything is animatable jow_go_for2_5 +++++++*/ /* introducing them here, because i know: steps in properties ( at frame timing ) * will cause unwanted responses of the softbody system (which does inter frame calculations ) * so first 'cure' would be: interpolate linear in time .. * Q: why do i write this? * A: because it happend once, that some eger coder 'streamlined' code to fail. * We DO linear interpolation for goals .. and i think we should do on animated properties as well */ /* animate sb->maxgoal, sb->mingoal */ static float _final_goal(Object *ob, BodyPoint *bp)/*jow_go_for2_5 */ { float f = -1999.99f; if (ob) { SoftBody *sb= ob->soft; /* is supposed to be there */ if (!(ob->softflag & OB_SB_GOAL)) return (0.0f); if (sb&&bp) { if (bp->goal < 0.0f) return (0.0f); f = sb->mingoal + bp->goal * fabsf(sb->maxgoal - sb->mingoal); f = pow(f, 4.0f); return (f); } } printf("_final_goal failed! sb or bp ==NULL\n"); return f; /*using crude but spot able values some times helps debuggin */ } static float _final_mass(Object *ob, BodyPoint *bp) { if (ob) { SoftBody *sb= ob->soft; /* is supposed to be there */ if (sb&&bp) { return(bp->mass*sb->nodemass); } } printf("_final_mass failed! sb or bp ==NULL\n"); return 1.0f; } /* helper functions for everything is animateble jow_go_for2_5 ------*/ /*+++ collider caching and dicing +++*/ /******************** for each target object/face the axis aligned bounding box (AABB) is stored faces parallel to global axes so only simple "value" in [min, max] ckecks are used float operations still */ /* just an ID here to reduce the prob for killing objects ** ob->sumohandle points to we should not kill :) */ static const int CCD_SAVETY = 190561; typedef struct ccdf_minmax { float minx, miny, minz, maxx, maxy, maxz; } ccdf_minmax; typedef struct ccd_Mesh { int totvert, totface; MVert *mvert; MVert *mprevvert; MFace *mface; int savety; ccdf_minmax *mima; /* Axis Aligned Bounding Box AABB */ float bbmin[3]; float bbmax[3]; } ccd_Mesh; static ccd_Mesh *ccd_mesh_make(Object *ob) { CollisionModifierData *cmd; ccd_Mesh *pccd_M = NULL; ccdf_minmax *mima = NULL; MFace *mface=NULL; float v[3], hull; int i; cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision); /* first some paranoia checks */ if (!cmd) return NULL; if (!cmd->numverts || !cmd->numfaces) return NULL; pccd_M = MEM_mallocN(sizeof(ccd_Mesh), "ccd_Mesh"); pccd_M->totvert = cmd->numverts; pccd_M->totface = cmd->numfaces; pccd_M->savety = CCD_SAVETY; pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f; pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f; pccd_M->mprevvert=NULL; /* blow it up with forcefield ranges */ hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft); /* alloc and copy verts*/ pccd_M->mvert = MEM_dupallocN(cmd->xnew); /* note that xnew coords are already in global space, */ /* determine the ortho BB */ for (i=0; i < pccd_M->totvert; i++) { /* evaluate limits */ copy_v3_v3(v, pccd_M->mvert[i].co); pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull); pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull); pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull); pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull); pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull); pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull); } /* alloc and copy faces*/ pccd_M->mface = MEM_dupallocN(cmd->mfaces); /* OBBs for idea1 */ pccd_M->mima = MEM_mallocN(sizeof(ccdf_minmax)*pccd_M->totface, "ccd_Mesh_Faces_mima"); mima = pccd_M->mima; mface = pccd_M->mface; /* anyhoo we need to walk the list of faces and find OBB they live in */ for (i=0; i < pccd_M->totface; i++) { mima->minx=mima->miny=mima->minz=1e30f; mima->maxx=mima->maxy=mima->maxz=-1e30f; copy_v3_v3(v, pccd_M->mvert[mface->v1].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mvert[mface->v2].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mvert[mface->v3].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); if (mface->v4) { copy_v3_v3(v, pccd_M->mvert[mface->v4].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); } mima++; mface++; } return pccd_M; } static void ccd_mesh_update(Object *ob, ccd_Mesh *pccd_M) { CollisionModifierData *cmd; ccdf_minmax *mima = NULL; MFace *mface=NULL; float v[3], hull; int i; cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision); /* first some paranoia checks */ if (!cmd) return; if (!cmd->numverts || !cmd->numfaces) return; if ((pccd_M->totvert != cmd->numverts) || (pccd_M->totface != cmd->numfaces)) { return; } pccd_M->bbmin[0]=pccd_M->bbmin[1]=pccd_M->bbmin[2]=1e30f; pccd_M->bbmax[0]=pccd_M->bbmax[1]=pccd_M->bbmax[2]=-1e30f; /* blow it up with forcefield ranges */ hull = max_ff(ob->pd->pdef_sbift, ob->pd->pdef_sboft); /* rotate current to previous */ if (pccd_M->mprevvert) MEM_freeN(pccd_M->mprevvert); pccd_M->mprevvert = pccd_M->mvert; /* alloc and copy verts*/ pccd_M->mvert = MEM_dupallocN(cmd->xnew); /* note that xnew coords are already in global space, */ /* determine the ortho BB */ for (i=0; i < pccd_M->totvert; i++) { /* evaluate limits */ copy_v3_v3(v, pccd_M->mvert[i].co); pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull); pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull); pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull); pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull); pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull); pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull); /* evaluate limits */ copy_v3_v3(v, pccd_M->mprevvert[i].co); pccd_M->bbmin[0] = min_ff(pccd_M->bbmin[0], v[0] - hull); pccd_M->bbmin[1] = min_ff(pccd_M->bbmin[1], v[1] - hull); pccd_M->bbmin[2] = min_ff(pccd_M->bbmin[2], v[2] - hull); pccd_M->bbmax[0] = max_ff(pccd_M->bbmax[0], v[0] + hull); pccd_M->bbmax[1] = max_ff(pccd_M->bbmax[1], v[1] + hull); pccd_M->bbmax[2] = max_ff(pccd_M->bbmax[2], v[2] + hull); } mima = pccd_M->mima; mface = pccd_M->mface; /* anyhoo we need to walk the list of faces and find OBB they live in */ for (i=0; i < pccd_M->totface; i++) { mima->minx=mima->miny=mima->minz=1e30f; mima->maxx=mima->maxy=mima->maxz=-1e30f; copy_v3_v3(v, pccd_M->mvert[mface->v1].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mvert[mface->v2].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mvert[mface->v3].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); if (mface->v4) { copy_v3_v3(v, pccd_M->mvert[mface->v4].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); } copy_v3_v3(v, pccd_M->mprevvert[mface->v1].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mprevvert[mface->v2].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); copy_v3_v3(v, pccd_M->mprevvert[mface->v3].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); if (mface->v4) { copy_v3_v3(v, pccd_M->mprevvert[mface->v4].co); mima->minx = min_ff(mima->minx, v[0] - hull); mima->miny = min_ff(mima->miny, v[1] - hull); mima->minz = min_ff(mima->minz, v[2] - hull); mima->maxx = max_ff(mima->maxx, v[0] + hull); mima->maxy = max_ff(mima->maxy, v[1] + hull); mima->maxz = max_ff(mima->maxz, v[2] + hull); } mima++; mface++; } return; } static void ccd_mesh_free(ccd_Mesh *ccdm) { if (ccdm && (ccdm->savety == CCD_SAVETY )) { /*make sure we're not nuking objects we don't know*/ MEM_freeN(ccdm->mface); MEM_freeN(ccdm->mvert); if (ccdm->mprevvert) MEM_freeN(ccdm->mprevvert); MEM_freeN(ccdm->mima); MEM_freeN(ccdm); ccdm = NULL; } } static void ccd_build_deflector_hash(Scene *scene, Object *vertexowner, GHash *hash) { Base *base= scene->base.first; Object *ob; if (!hash) return; while (base) { /*Only proceed for mesh object in same layer */ if (base->object->type==OB_MESH && (base->lay & vertexowner->lay)) { ob= base->object; if ((vertexowner) && (ob == vertexowner)) { /* if vertexowner is given we don't want to check collision with owner object */ base = base->next; continue; } /*+++ only with deflecting set */ if (ob->pd && ob->pd->deflect && BLI_ghash_lookup(hash, ob) == NULL) { ccd_Mesh *ccdmesh = ccd_mesh_make(ob); BLI_ghash_insert(hash, ob, ccdmesh); }/*--- only with deflecting set */ }/* mesh && layer*/ base = base->next; } /* while (base) */ } static void ccd_update_deflector_hash(Scene *scene, Object *vertexowner, GHash *hash) { Base *base= scene->base.first; Object *ob; if ((!hash) || (!vertexowner)) return; while (base) { /*Only proceed for mesh object in same layer */ if (base->object->type==OB_MESH && (base->lay & vertexowner->lay)) { ob= base->object; if (ob == vertexowner) { /* if vertexowner is given we don't want to check collision with owner object */ base = base->next; continue; } /*+++ only with deflecting set */ if (ob->pd && ob->pd->deflect) { ccd_Mesh *ccdmesh = BLI_ghash_lookup(hash, ob); if (ccdmesh) ccd_mesh_update(ob, ccdmesh); }/*--- only with deflecting set */ }/* mesh && layer*/ base = base->next; } /* while (base) */ } /*--- collider caching and dicing ---*/ static int count_mesh_quads(Mesh *me) { int a, result = 0; MFace *mface= me->mface; if (mface) { for (a=me->totface; a>0; a--, mface++) { if (mface->v4) result++; } } return result; } static void add_mesh_quad_diag_springs(Object *ob) { Mesh *me= ob->data; MFace *mface= me->mface; /*BodyPoint *bp;*/ /*UNUSED*/ BodySpring *bs, *bs_new; int a; if (ob->soft) { int nofquads; //float s_shear = ob->soft->shearstiff*ob->soft->shearstiff; nofquads = count_mesh_quads(me); if (nofquads) { /* resize spring-array to hold additional quad springs */ bs_new= MEM_callocN((ob->soft->totspring + nofquads *2 )*sizeof(BodySpring), "bodyspring"); memcpy(bs_new, ob->soft->bspring, (ob->soft->totspring )*sizeof(BodySpring)); if (ob->soft->bspring) MEM_freeN(ob->soft->bspring); /* do this before reassigning the pointer or have a 1st class memory leak */ ob->soft->bspring = bs_new; /* fill the tail */ a = 0; bs = bs_new+ob->soft->totspring; /*bp= ob->soft->bpoint; */ /*UNUSED*/ if (mface ) { for (a=me->totface; a>0; a--, mface++) { if (mface->v4) { bs->v1= mface->v1; bs->v2= mface->v3; bs->springtype = SB_STIFFQUAD; bs++; bs->v1= mface->v2; bs->v2= mface->v4; bs->springtype = SB_STIFFQUAD; bs++; } } } /* now we can announce new springs */ ob->soft->totspring += nofquads *2; } } } static void add_2nd_order_roller(Object *ob, float UNUSED(stiffness), int *counter, int addsprings) { /*assume we have a softbody*/ SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp, *bpo; BodySpring *bs, *bs2, *bs3= NULL; int a, b, c, notthis= 0, v0; if (!sb->bspring) {return;} /* we are 2nd order here so 1rst should have been build :) */ /* first run counting second run adding */ *counter = 0; if (addsprings) bs3 = ob->soft->bspring+ob->soft->totspring; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { /*scan for neighborhood*/ bpo = NULL; v0 = (sb->totpoint-a); for (b=bp->nofsprings;b>0;b--) { bs = sb->bspring + bp->springs[b-1]; /*nasty thing here that springs have two ends so here we have to make sure we examine the other */ if (v0 == bs->v1) { bpo = sb->bpoint+bs->v2; notthis = bs->v2; } else { if (v0 == bs->v2) { bpo = sb->bpoint+bs->v1; notthis = bs->v1; } else {printf("oops we should not get here - add_2nd_order_springs");} } if (bpo) {/* so now we have a 2nd order humpdidump */ for (c=bpo->nofsprings;c>0;c--) { bs2 = sb->bspring + bpo->springs[c-1]; if ((bs2->v1 != notthis) && (bs2->v1 > v0)) { (*counter)++;/*hit */ if (addsprings) { bs3->v1= v0; bs3->v2= bs2->v1; bs3->springtype = SB_BEND; bs3++; } } if ((bs2->v2 !=notthis)&&(bs2->v2 > v0)) { (*counter)++; /* hit */ if (addsprings) { bs3->v1= v0; bs3->v2= bs2->v2; bs3->springtype = SB_BEND; bs3++; } } } } } /*scan for neighborhood done*/ } } static void add_2nd_order_springs(Object *ob, float stiffness) { int counter = 0; BodySpring *bs_new; stiffness *=stiffness; add_2nd_order_roller(ob, stiffness, &counter, 0); /* counting */ if (counter) { /* resize spring-array to hold additional springs */ bs_new= MEM_callocN((ob->soft->totspring + counter )*sizeof(BodySpring), "bodyspring"); memcpy(bs_new, ob->soft->bspring, (ob->soft->totspring )*sizeof(BodySpring)); if (ob->soft->bspring) MEM_freeN(ob->soft->bspring); ob->soft->bspring = bs_new; add_2nd_order_roller(ob, stiffness, &counter, 1); /* adding */ ob->soft->totspring += counter; } } static void add_bp_springlist(BodyPoint *bp, int springID) { int *newlist; if (bp->springs == NULL) { bp->springs = MEM_callocN(sizeof(int), "bpsprings"); bp->springs[0] = springID; bp->nofsprings = 1; } else { bp->nofsprings++; newlist = MEM_callocN(bp->nofsprings * sizeof(int), "bpsprings"); memcpy(newlist, bp->springs, (bp->nofsprings-1)* sizeof(int)); MEM_freeN(bp->springs); bp->springs = newlist; bp->springs[bp->nofsprings-1] = springID; } } /* do this once when sb is build it is O(N^2) so scanning for springs every iteration is too expensive */ static void build_bps_springlist(Object *ob) { SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; BodySpring *bs; int a, b; if (sb==NULL) return; /* paranoya check */ for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { /* throw away old list */ if (bp->springs) { MEM_freeN(bp->springs); bp->springs=NULL; } /* scan for attached inner springs */ for (b=sb->totspring, bs= sb->bspring; b>0; b--, bs++) { if (( (sb->totpoint-a) == bs->v1) ) { add_bp_springlist(bp, sb->totspring -b); } if (( (sb->totpoint-a) == bs->v2) ) { add_bp_springlist(bp, sb->totspring -b); } }/*for springs*/ }/*for bp*/ } static void calculate_collision_balls(Object *ob) { SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; BodySpring *bs; int a, b, akku_count; float min, max, akku; if (sb==NULL) return; /* paranoya check */ for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { bp->colball=0; akku =0.0f; akku_count=0; min = 1e22f; max = -1e22f; /* first estimation based on attached */ for (b=bp->nofsprings;b>0;b--) { bs = sb->bspring + bp->springs[b-1]; if (bs->springtype == SB_EDGE) { akku += bs->len; akku_count++; min = min_ff(bs->len, min); max = max_ff(bs->len, max); } } if (akku_count > 0) { if (sb->sbc_mode == SBC_MODE_MANUAL) { bp->colball=sb->colball; } if (sb->sbc_mode == SBC_MODE_AVG) { bp->colball = akku/(float)akku_count*sb->colball; } if (sb->sbc_mode == SBC_MODE_MIN) { bp->colball=min*sb->colball; } if (sb->sbc_mode == SBC_MODE_MAX) { bp->colball=max*sb->colball; } if (sb->sbc_mode == SBC_MODE_AVGMINMAX) { bp->colball = (min + max)/2.0f*sb->colball; } } else bp->colball=0; }/*for bp*/ } /* creates new softbody if didn't exist yet, makes new points and springs arrays */ static void renew_softbody(Scene *scene, Object *ob, int totpoint, int totspring) { SoftBody *sb; int i; short softflag; if (ob->soft==NULL) ob->soft= sbNew(scene); else free_softbody_intern(ob->soft); sb= ob->soft; softflag=ob->softflag; if (totpoint) { sb->totpoint= totpoint; sb->totspring= totspring; sb->bpoint= MEM_mallocN(totpoint*sizeof(BodyPoint), "bodypoint"); if (totspring) sb->bspring= MEM_mallocN(totspring*sizeof(BodySpring), "bodyspring"); /* initialize BodyPoint array */ for (i=0; ibpoint[i]; /* hum as far as i see this is overridden by _final_goal() now jow_go_for2_5 */ /* sadly breaks compatibility with older versions */ /* but makes goals behave the same for meshes, lattices and curves */ if (softflag & OB_SB_GOAL) { bp->goal= sb->defgoal; } else { bp->goal= 0.0f; /* so this will definily be below SOFTGOALSNAP */ } bp->nofsprings= 0; bp->springs= NULL; bp->choke = 0.0f; bp->choke2 = 0.0f; bp->frozen = 1.0f; bp->colball = 0.0f; bp->loc_flag = 0; bp->springweight = 1.0f; bp->mass = 1.0f; } } } static void free_softbody_baked(SoftBody *sb) { SBVertex *key; int k; for (k=0; ktotkey; k++) { key= *(sb->keys + k); if (key) MEM_freeN(key); } if (sb->keys) MEM_freeN(sb->keys); sb->keys= NULL; sb->totkey= 0; } static void free_scratch(SoftBody *sb) { if (sb->scratch) { /* todo make sure everything is cleaned up nicly */ if (sb->scratch->colliderhash) { BLI_ghash_free(sb->scratch->colliderhash, NULL, (GHashValFreeFP) ccd_mesh_free); /*this hoepfully will free all caches*/ sb->scratch->colliderhash = NULL; } if (sb->scratch->bodyface) { MEM_freeN(sb->scratch->bodyface); } if (sb->scratch->Ref.ivert) { MEM_freeN(sb->scratch->Ref.ivert); } MEM_freeN(sb->scratch); sb->scratch = NULL; } } /* only frees internal data */ static void free_softbody_intern(SoftBody *sb) { if (sb) { int a; BodyPoint *bp; if (sb->bpoint) { for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { /* free spring list */ if (bp->springs != NULL) { MEM_freeN(bp->springs); } } MEM_freeN(sb->bpoint); } if (sb->bspring) MEM_freeN(sb->bspring); sb->totpoint= sb->totspring= 0; sb->bpoint= NULL; sb->bspring= NULL; free_scratch(sb); free_softbody_baked(sb); } } /* ************ dynamics ********** */ /* the most general (micro physics correct) way to do collision ** (only needs the current particle position) ** ** it actually checks if the particle intrudes a short range force field generated ** by the faces of the target object and returns a force to drive the particel out ** the strenght of the field grows exponetially if the particle is on the 'wrong' side of the face ** 'wrong' side : projection to the face normal is negative (all referred to a vertex in the face) ** ** flaw of this: 'fast' particles as well as 'fast' colliding faces ** give a 'tunnel' effect such that the particle passes through the force field ** without ever 'seeing' it ** this is fully compliant to heisenberg: h >= fuzzy(location) * fuzzy(time) ** besides our h is way larger than in QM because forces propagate way slower here ** we have to deal with fuzzy(time) in the range of 1/25 seconds (typical frame rate) ** yup collision targets are not known here any better ** and 1/25 second is looong compared to real collision events ** Q: why not use 'simple' collision here like bouncing back a particle ** --> reverting is velocity on the face normal ** A: because our particles are not alone here ** and need to tell their neighbors exactly what happens via spring forces ** unless sbObjectStep( .. ) is called on sub frame timing level ** BTW that also questions the use of a 'implicit' solvers on softbodies ** since that would only valid for 'slow' moving collision targets and dito particles */ /* aye this belongs to arith.c */ static void Vec3PlusStVec(float *v, float s, float *v1) { v[0] += s*v1[0]; v[1] += s*v1[1]; v[2] += s*v1[2]; } /* +++ dependency information functions*/ static int are_there_deflectors(Scene *scene, unsigned int layer) { Base *base; for (base = scene->base.first; base; base= base->next) { if ( (base->lay & layer) && base->object->pd) { if (base->object->pd->deflect) return 1; } } return 0; } static int query_external_colliders(Scene *scene, Object *me) { return(are_there_deflectors(scene, me->lay)); } /* --- dependency information functions*/ /* +++ the aabb "force" section*/ static int sb_detect_aabb_collisionCached(float UNUSED(force[3]), unsigned int UNUSED(par_layer), struct Object *vertexowner, float UNUSED(time)) { Object *ob; SoftBody *sb=vertexowner->soft; GHash *hash; GHashIterator *ihash; float aabbmin[3], aabbmax[3]; int deflected=0; #if 0 int a; #endif if ((sb == NULL) || (sb->scratch ==NULL)) return 0; copy_v3_v3(aabbmin, sb->scratch->aabbmin); copy_v3_v3(aabbmax, sb->scratch->aabbmax); hash = vertexowner->soft->scratch->colliderhash; ihash = BLI_ghashIterator_new(hash); while (!BLI_ghashIterator_done(ihash)) { ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); ob = BLI_ghashIterator_getKey (ihash); { /* only with deflecting set */ if (ob->pd && ob->pd->deflect) { #if 0 /* UNUSED */ MFace *mface= NULL; MVert *mvert= NULL; MVert *mprevvert= NULL; ccdf_minmax *mima= NULL; #endif if (ccdm) { #if 0 /* UNUSED */ mface= ccdm->mface; mvert= ccdm->mvert; mprevvert= ccdm->mprevvert; mima= ccdm->mima; a = ccdm->totface; #endif if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) || (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) || (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2]) ) { /* boxes don't intersect */ BLI_ghashIterator_step(ihash); continue; } /* so now we have the 2 boxes overlapping */ /* forces actually not used */ deflected = 2; } else { /*aye that should be cached*/ printf("missing cache error\n"); BLI_ghashIterator_step(ihash); continue; } } /* if (ob->pd && ob->pd->deflect) */ BLI_ghashIterator_step(ihash); } } /* while () */ BLI_ghashIterator_free(ihash); return deflected; } /* --- the aabb section*/ /* +++ the face external section*/ static int sb_detect_face_pointCached(float face_v1[3], float face_v2[3], float face_v3[3], float *damp, float force[3], unsigned int UNUSED(par_layer), struct Object *vertexowner, float time) { Object *ob; GHash *hash; GHashIterator *ihash; float nv1[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3]; float facedist, outerfacethickness, tune = 10.f; int a, deflected=0; aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]); aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]); aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]); aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]); aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]); aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]); /* calculate face normal once again SIGH */ sub_v3_v3v3(edge1, face_v1, face_v2); sub_v3_v3v3(edge2, face_v3, face_v2); cross_v3_v3v3(d_nvect, edge2, edge1); normalize_v3(d_nvect); hash = vertexowner->soft->scratch->colliderhash; ihash = BLI_ghashIterator_new(hash); while (!BLI_ghashIterator_done(ihash)) { ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); ob = BLI_ghashIterator_getKey (ihash); { /* only with deflecting set */ if (ob->pd && ob->pd->deflect) { MVert *mvert= NULL; MVert *mprevvert= NULL; if (ccdm) { mvert= ccdm->mvert; a = ccdm->totvert; mprevvert= ccdm->mprevvert; outerfacethickness = ob->pd->pdef_sboft; if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) || (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) || (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2]) ) { /* boxes don't intersect */ BLI_ghashIterator_step(ihash); continue; } } else { /*aye that should be cached*/ printf("missing cache error\n"); BLI_ghashIterator_step(ihash); continue; } /* use mesh*/ if (mvert) { while (a) { copy_v3_v3(nv1, mvert[a-1].co); if (mprevvert) { mul_v3_fl(nv1, time); Vec3PlusStVec(nv1, (1.0f-time), mprevvert[a-1].co); } /* origin to face_v2*/ sub_v3_v3(nv1, face_v2); facedist = dot_v3v3(nv1, d_nvect); if (ABS(facedist) 0) { df = (outerfacethickness-facedist)/outerfacethickness; } else { df = (outerfacethickness+facedist)/outerfacethickness; } *damp=df*tune*ob->pd->pdef_sbdamp; df = 0.01f * expf(-100.0f * df); Vec3PlusStVec(force, -df, d_nvect); deflected = 3; } } a--; }/* while (a)*/ } /* if (mvert) */ } /* if (ob->pd && ob->pd->deflect) */ BLI_ghashIterator_step(ihash); } } /* while () */ BLI_ghashIterator_free(ihash); return deflected; } static int sb_detect_face_collisionCached(float face_v1[3], float face_v2[3], float face_v3[3], float *damp, float force[3], unsigned int UNUSED(par_layer), struct Object *vertexowner, float time) { Object *ob; GHash *hash; GHashIterator *ihash; float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3]; float t, tune = 10.0f; int a, deflected=0; aabbmin[0] = min_fff(face_v1[0], face_v2[0], face_v3[0]); aabbmin[1] = min_fff(face_v1[1], face_v2[1], face_v3[1]); aabbmin[2] = min_fff(face_v1[2], face_v2[2], face_v3[2]); aabbmax[0] = max_fff(face_v1[0], face_v2[0], face_v3[0]); aabbmax[1] = max_fff(face_v1[1], face_v2[1], face_v3[1]); aabbmax[2] = max_fff(face_v1[2], face_v2[2], face_v3[2]); hash = vertexowner->soft->scratch->colliderhash; ihash = BLI_ghashIterator_new(hash); while (!BLI_ghashIterator_done(ihash)) { ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); ob = BLI_ghashIterator_getKey (ihash); { /* only with deflecting set */ if (ob->pd && ob->pd->deflect) { MFace *mface= NULL; MVert *mvert= NULL; MVert *mprevvert= NULL; ccdf_minmax *mima= NULL; if (ccdm) { mface= ccdm->mface; mvert= ccdm->mvert; mprevvert= ccdm->mprevvert; mima= ccdm->mima; a = ccdm->totface; if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) || (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) || (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2]) ) { /* boxes don't intersect */ BLI_ghashIterator_step(ihash); continue; } } else { /*aye that should be cached*/ printf("missing cache error\n"); BLI_ghashIterator_step(ihash); continue; } /* use mesh*/ while (a--) { if ( (aabbmax[0] < mima->minx) || (aabbmin[0] > mima->maxx) || (aabbmax[1] < mima->miny) || (aabbmin[1] > mima->maxy) || (aabbmax[2] < mima->minz) || (aabbmin[2] > mima->maxz) ) { mface++; mima++; continue; } if (mvert) { copy_v3_v3(nv1, mvert[mface->v1].co); copy_v3_v3(nv2, mvert[mface->v2].co); copy_v3_v3(nv3, mvert[mface->v3].co); if (mface->v4) { copy_v3_v3(nv4, mvert[mface->v4].co); } if (mprevvert) { mul_v3_fl(nv1, time); Vec3PlusStVec(nv1, (1.0f-time), mprevvert[mface->v1].co); mul_v3_fl(nv2, time); Vec3PlusStVec(nv2, (1.0f-time), mprevvert[mface->v2].co); mul_v3_fl(nv3, time); Vec3PlusStVec(nv3, (1.0f-time), mprevvert[mface->v3].co); if (mface->v4) { mul_v3_fl(nv4, time); Vec3PlusStVec(nv4, (1.0f-time), mprevvert[mface->v4].co); } } } /* switch origin to be nv2*/ sub_v3_v3v3(edge1, nv1, nv2); sub_v3_v3v3(edge2, nv3, nv2); cross_v3_v3v3(d_nvect, edge2, edge1); normalize_v3(d_nvect); if ( isect_line_tri_v3(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) || isect_line_tri_v3(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) || isect_line_tri_v3(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL) ) { Vec3PlusStVec(force, -0.5f, d_nvect); *damp=tune*ob->pd->pdef_sbdamp; deflected = 2; } if (mface->v4) { /* quad */ /* switch origin to be nv4 */ sub_v3_v3v3(edge1, nv3, nv4); sub_v3_v3v3(edge2, nv1, nv4); cross_v3_v3v3(d_nvect, edge2, edge1); normalize_v3(d_nvect); if ( /* isect_line_tri_v3(nv1, nv3, face_v1, face_v2, face_v3, &t, NULL) || * we did that edge already */ isect_line_tri_v3(nv3, nv4, face_v1, face_v2, face_v3, &t, NULL) || isect_line_tri_v3(nv4, nv1, face_v1, face_v2, face_v3, &t, NULL) ) { Vec3PlusStVec(force, -0.5f, d_nvect); *damp=tune*ob->pd->pdef_sbdamp; deflected = 2; } } mface++; mima++; }/* while a */ } /* if (ob->pd && ob->pd->deflect) */ BLI_ghashIterator_step(ihash); } } /* while () */ BLI_ghashIterator_free(ihash); return deflected; } static void scan_for_ext_face_forces(Object *ob, float timenow) { SoftBody *sb = ob->soft; BodyFace *bf; int a; float damp=0.0f, choke=1.0f; float tune = -10.0f; float feedback[3]; if (sb && sb->scratch->totface) { bf = sb->scratch->bodyface; for (a=0; ascratch->totface; a++, bf++) { bf->ext_force[0]=bf->ext_force[1]=bf->ext_force[2]=0.0f; /*+++edges intruding*/ bf->flag &= ~BFF_INTERSECT; feedback[0]=feedback[1]=feedback[2]=0.0f; if (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, &damp, feedback, ob->lay, ob, timenow)) { Vec3PlusStVec(sb->bpoint[bf->v1].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v2].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v3].force, tune, feedback); // Vec3PlusStVec(bf->ext_force, tune, feedback); bf->flag |= BFF_INTERSECT; choke = min_ff(max_ff(damp, choke), 1.0f); } feedback[0]=feedback[1]=feedback[2]=0.0f; if ((bf->v4) && (sb_detect_face_collisionCached(sb->bpoint[bf->v1].pos, sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, &damp, feedback, ob->lay, ob, timenow))) { Vec3PlusStVec(sb->bpoint[bf->v1].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v3].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v4].force, tune, feedback); // Vec3PlusStVec(bf->ext_force, tune, feedback); bf->flag |= BFF_INTERSECT; choke = min_ff(max_ff(damp, choke), 1.0f); } /*---edges intruding*/ /*+++ close vertices*/ if (( bf->flag & BFF_INTERSECT)==0) { bf->flag &= ~BFF_CLOSEVERT; tune = -1.0f; feedback[0]=feedback[1]=feedback[2]=0.0f; if (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos, &damp, feedback, ob->lay, ob, timenow)) { Vec3PlusStVec(sb->bpoint[bf->v1].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v2].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v3].force, tune, feedback); // Vec3PlusStVec(bf->ext_force, tune, feedback); bf->flag |= BFF_CLOSEVERT; choke = min_ff(max_ff(damp, choke), 1.0f); } feedback[0]=feedback[1]=feedback[2]=0.0f; if ((bf->v4) && (sb_detect_face_pointCached(sb->bpoint[bf->v1].pos, sb->bpoint[bf->v3].pos, sb->bpoint[bf->v4].pos, &damp, feedback, ob->lay, ob, timenow))) { Vec3PlusStVec(sb->bpoint[bf->v1].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v3].force, tune, feedback); Vec3PlusStVec(sb->bpoint[bf->v4].force, tune, feedback); // Vec3PlusStVec(bf->ext_force, tune, feedback); bf->flag |= BFF_CLOSEVERT; choke = min_ff(max_ff(damp, choke), 1.0f); } } /*--- close vertices*/ } bf = sb->scratch->bodyface; for (a=0; ascratch->totface; a++, bf++) { if (( bf->flag & BFF_INTERSECT) || ( bf->flag & BFF_CLOSEVERT)) { sb->bpoint[bf->v1].choke2 = max_ff(sb->bpoint[bf->v1].choke2, choke); sb->bpoint[bf->v2].choke2 = max_ff(sb->bpoint[bf->v2].choke2, choke); sb->bpoint[bf->v3].choke2 = max_ff(sb->bpoint[bf->v3].choke2, choke); if (bf->v4) { sb->bpoint[bf->v2].choke2 = max_ff(sb->bpoint[bf->v2].choke2, choke); } } } } } /* --- the face external section*/ /* +++ the spring external section*/ static int sb_detect_edge_collisionCached(float edge_v1[3], float edge_v2[3], float *damp, float force[3], unsigned int UNUSED(par_layer), struct Object *vertexowner, float time) { Object *ob; GHash *hash; GHashIterator *ihash; float nv1[3], nv2[3], nv3[3], nv4[3], edge1[3], edge2[3], d_nvect[3], aabbmin[3], aabbmax[3]; float t, el; int a, deflected=0; minmax_v3v3_v3(aabbmin, aabbmax, edge_v1); minmax_v3v3_v3(aabbmin, aabbmax, edge_v2); el = len_v3v3(edge_v1, edge_v2); hash = vertexowner->soft->scratch->colliderhash; ihash = BLI_ghashIterator_new(hash); while (!BLI_ghashIterator_done(ihash)) { ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); ob = BLI_ghashIterator_getKey (ihash); { /* only with deflecting set */ if (ob->pd && ob->pd->deflect) { MFace *mface= NULL; MVert *mvert= NULL; MVert *mprevvert= NULL; ccdf_minmax *mima= NULL; if (ccdm) { mface= ccdm->mface; mvert= ccdm->mvert; mprevvert= ccdm->mprevvert; mima= ccdm->mima; a = ccdm->totface; if ((aabbmax[0] < ccdm->bbmin[0]) || (aabbmax[1] < ccdm->bbmin[1]) || (aabbmax[2] < ccdm->bbmin[2]) || (aabbmin[0] > ccdm->bbmax[0]) || (aabbmin[1] > ccdm->bbmax[1]) || (aabbmin[2] > ccdm->bbmax[2]) ) { /* boxes don't intersect */ BLI_ghashIterator_step(ihash); continue; } } else { /*aye that should be cached*/ printf("missing cache error\n"); BLI_ghashIterator_step(ihash); continue; } /* use mesh*/ while (a--) { if ( (aabbmax[0] < mima->minx) || (aabbmin[0] > mima->maxx) || (aabbmax[1] < mima->miny) || (aabbmin[1] > mima->maxy) || (aabbmax[2] < mima->minz) || (aabbmin[2] > mima->maxz) ) { mface++; mima++; continue; } if (mvert) { copy_v3_v3(nv1, mvert[mface->v1].co); copy_v3_v3(nv2, mvert[mface->v2].co); copy_v3_v3(nv3, mvert[mface->v3].co); if (mface->v4) { copy_v3_v3(nv4, mvert[mface->v4].co); } if (mprevvert) { mul_v3_fl(nv1, time); Vec3PlusStVec(nv1, (1.0f-time), mprevvert[mface->v1].co); mul_v3_fl(nv2, time); Vec3PlusStVec(nv2, (1.0f-time), mprevvert[mface->v2].co); mul_v3_fl(nv3, time); Vec3PlusStVec(nv3, (1.0f-time), mprevvert[mface->v3].co); if (mface->v4) { mul_v3_fl(nv4, time); Vec3PlusStVec(nv4, (1.0f-time), mprevvert[mface->v4].co); } } } /* switch origin to be nv2*/ sub_v3_v3v3(edge1, nv1, nv2); sub_v3_v3v3(edge2, nv3, nv2); cross_v3_v3v3(d_nvect, edge2, edge1); normalize_v3(d_nvect); if ( isect_line_tri_v3(edge_v1, edge_v2, nv1, nv2, nv3, &t, NULL)) { float v1[3], v2[3]; float intrusiondepth, i1, i2; sub_v3_v3v3(v1, edge_v1, nv2); sub_v3_v3v3(v2, edge_v2, nv2); i1 = dot_v3v3(v1, d_nvect); i2 = dot_v3v3(v2, d_nvect); intrusiondepth = -min_ff(i1, i2) / el; Vec3PlusStVec(force, intrusiondepth, d_nvect); *damp=ob->pd->pdef_sbdamp; deflected = 2; } if (mface->v4) { /* quad */ /* switch origin to be nv4 */ sub_v3_v3v3(edge1, nv3, nv4); sub_v3_v3v3(edge2, nv1, nv4); cross_v3_v3v3(d_nvect, edge2, edge1); normalize_v3(d_nvect); if (isect_line_tri_v3( edge_v1, edge_v2, nv1, nv3, nv4, &t, NULL)) { float v1[3], v2[3]; float intrusiondepth, i1, i2; sub_v3_v3v3(v1, edge_v1, nv4); sub_v3_v3v3(v2, edge_v2, nv4); i1 = dot_v3v3(v1, d_nvect); i2 = dot_v3v3(v2, d_nvect); intrusiondepth = -min_ff(i1, i2) / el; Vec3PlusStVec(force, intrusiondepth, d_nvect); *damp=ob->pd->pdef_sbdamp; deflected = 2; } } mface++; mima++; }/* while a */ } /* if (ob->pd && ob->pd->deflect) */ BLI_ghashIterator_step(ihash); } } /* while () */ BLI_ghashIterator_free(ihash); return deflected; } static void _scan_for_ext_spring_forces(Scene *scene, Object *ob, float timenow, int ifirst, int ilast, struct ListBase *do_effector) { SoftBody *sb = ob->soft; int a; float damp; float feedback[3]; if (sb && sb->totspring) { for (a=ifirst; abspring[a]; bs->ext_force[0]=bs->ext_force[1]=bs->ext_force[2]=0.0f; feedback[0]=feedback[1]=feedback[2]=0.0f; bs->flag &= ~BSF_INTERSECT; if (bs->springtype == SB_EDGE) { /* +++ springs colliding */ if (ob->softflag & OB_SB_EDGECOLL) { if ( sb_detect_edge_collisionCached (sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos, &damp, feedback, ob->lay, ob, timenow)) { add_v3_v3(bs->ext_force, feedback); bs->flag |= BSF_INTERSECT; //bs->cf=damp; bs->cf=sb->choke*0.01f; } } /* ---- springs colliding */ /* +++ springs seeing wind ... n stuff depending on their orientation*/ /* note we don't use sb->mediafrict but use sb->aeroedge for magnitude of effect*/ if (sb->aeroedge) { float vel[3], sp[3], pr[3], force[3]; float f, windfactor = 0.25f; /*see if we have wind*/ if (do_effector) { EffectedPoint epoint; float speed[3] = {0.0f, 0.0f, 0.0f}; float pos[3]; mid_v3_v3v3(pos, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos); mid_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec); pd_point_from_soft(scene, pos, vel, -1, &epoint); pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); mul_v3_fl(speed, windfactor); add_v3_v3(vel, speed); } /* media in rest */ else { add_v3_v3v3(vel, sb->bpoint[bs->v1].vec, sb->bpoint[bs->v2].vec); } f = normalize_v3(vel); f = -0.0001f*f*f*sb->aeroedge; /* (todo) add a nice angle dependent function done for now BUT */ /* still there could be some nice drag/lift function, but who needs it */ sub_v3_v3v3(sp, sb->bpoint[bs->v1].pos, sb->bpoint[bs->v2].pos); project_v3_v3v3(pr, vel, sp); sub_v3_v3(vel, pr); normalize_v3(vel); if (ob->softflag & OB_SB_AERO_ANGLE) { normalize_v3(sp); Vec3PlusStVec(bs->ext_force, f * (1.0f - fabsf(dot_v3v3(vel, sp))), vel); } else { Vec3PlusStVec(bs->ext_force, f, vel); // to keep compatible with 2.45 release files } } /* --- springs seeing wind */ } } } } static void scan_for_ext_spring_forces(Scene *scene, Object *ob, float timenow) { SoftBody *sb = ob->soft; ListBase *do_effector = NULL; do_effector = pdInitEffectors(scene, ob, NULL, sb->effector_weights, true); _scan_for_ext_spring_forces(scene, ob, timenow, 0, sb->totspring, do_effector); pdEndEffectors(&do_effector); } static void *exec_scan_for_ext_spring_forces(void *data) { SB_thread_context *pctx = (SB_thread_context*)data; _scan_for_ext_spring_forces(pctx->scene, pctx->ob, pctx->timenow, pctx->ifirst, pctx->ilast, pctx->do_effector); return NULL; } static void sb_sfesf_threads_run(Scene *scene, struct Object *ob, float timenow, int totsprings, int *UNUSED(ptr_to_break_func(void))) { ListBase *do_effector = NULL; ListBase threads; SB_thread_context *sb_threads; int i, totthread, left, dec; int lowsprings =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ do_effector= pdInitEffectors(scene, ob, NULL, ob->soft->effector_weights, true); /* figure the number of threads while preventing pretty pointless threading overhead */ totthread= BKE_scene_num_threads(scene); /* what if we got zillions of CPUs running but less to spread*/ while ((totsprings/totthread < lowsprings) && (totthread > 1)) { totthread--; } sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBSpringsThread"); memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); left = totsprings; dec = totsprings/totthread +1; for (i=0; i0) { sb_threads[i].ifirst = left; } else sb_threads[i].ifirst = 0; sb_threads[i].do_effector = do_effector; sb_threads[i].do_deflector = false;// not used here sb_threads[i].fieldfactor = 0.0f;// not used here sb_threads[i].windfactor = 0.0f;// not used here sb_threads[i].nr= i; sb_threads[i].tot= totthread; } if (totthread > 1) { BLI_init_threads(&threads, exec_scan_for_ext_spring_forces, totthread); for (i=0; isoft->scratch->colliderhash; ihash = BLI_ghashIterator_new(hash); outerforceaccu[0]=outerforceaccu[1]=outerforceaccu[2]=0.0f; innerforceaccu[0]=innerforceaccu[1]=innerforceaccu[2]=0.0f; /* go */ while (!BLI_ghashIterator_done(ihash)) { ccd_Mesh *ccdm = BLI_ghashIterator_getValue (ihash); ob = BLI_ghashIterator_getKey (ihash); { /* only with deflecting set */ if (ob->pd && ob->pd->deflect) { MFace *mface= NULL; MVert *mvert= NULL; MVert *mprevvert= NULL; ccdf_minmax *mima= NULL; if (ccdm) { mface= ccdm->mface; mvert= ccdm->mvert; mprevvert= ccdm->mprevvert; mima= ccdm->mima; a = ccdm->totface; minx = ccdm->bbmin[0]; miny = ccdm->bbmin[1]; minz = ccdm->bbmin[2]; maxx = ccdm->bbmax[0]; maxy = ccdm->bbmax[1]; maxz = ccdm->bbmax[2]; if ((opco[0] < minx) || (opco[1] < miny) || (opco[2] < minz) || (opco[0] > maxx) || (opco[1] > maxy) || (opco[2] > maxz) ) { /* outside the padded boundbox --> collision object is too far away */ BLI_ghashIterator_step(ihash); continue; } } else { /*aye that should be cached*/ printf("missing cache error\n"); BLI_ghashIterator_step(ihash); continue; } /* do object level stuff */ /* need to have user control for that since it depends on model scale */ innerfacethickness = -ob->pd->pdef_sbift; outerfacethickness = ob->pd->pdef_sboft; fa = (ff*outerfacethickness-outerfacethickness); fa *= fa; fa = 1.0f/fa; avel[0]=avel[1]=avel[2]=0.0f; /* use mesh*/ while (a--) { if ( (opco[0] < mima->minx) || (opco[0] > mima->maxx) || (opco[1] < mima->miny) || (opco[1] > mima->maxy) || (opco[2] < mima->minz) || (opco[2] > mima->maxz) ) { mface++; mima++; continue; } if (mvert) { copy_v3_v3(nv1, mvert[mface->v1].co); copy_v3_v3(nv2, mvert[mface->v2].co); copy_v3_v3(nv3, mvert[mface->v3].co); if (mface->v4) { copy_v3_v3(nv4, mvert[mface->v4].co); } if (mprevvert) { /* grab the average speed of the collider vertices before we spoil nvX humm could be done once a SB steps but then we' need to store that too since the AABB reduced propabitlty to get here drasticallly it might be a nice tradeof CPU <--> memory */ sub_v3_v3v3(vv1, nv1, mprevvert[mface->v1].co); sub_v3_v3v3(vv2, nv2, mprevvert[mface->v2].co); sub_v3_v3v3(vv3, nv3, mprevvert[mface->v3].co); if (mface->v4) { sub_v3_v3v3(vv4, nv4, mprevvert[mface->v4].co); } mul_v3_fl(nv1, time); Vec3PlusStVec(nv1, (1.0f-time), mprevvert[mface->v1].co); mul_v3_fl(nv2, time); Vec3PlusStVec(nv2, (1.0f-time), mprevvert[mface->v2].co); mul_v3_fl(nv3, time); Vec3PlusStVec(nv3, (1.0f-time), mprevvert[mface->v3].co); if (mface->v4) { mul_v3_fl(nv4, time); Vec3PlusStVec(nv4, (1.0f-time), mprevvert[mface->v4].co); } } } /* switch origin to be nv2*/ sub_v3_v3v3(edge1, nv1, nv2); sub_v3_v3v3(edge2, nv3, nv2); sub_v3_v3v3(dv1, opco, nv2); /* abuse dv1 to have vertex in question at *origin* of triangle */ cross_v3_v3v3(d_nvect, edge2, edge1); /* n_mag = */ /* UNUSED */ normalize_v3(d_nvect); facedist = dot_v3v3(dv1, d_nvect); // so rules are // if ((facedist > innerfacethickness) && (facedist < outerfacethickness)) { if (isect_point_tri_prism_v3(opco, nv1, nv2, nv3) ) { force_mag_norm =(float)exp(-ee*facedist); if (facedist > outerfacethickness*ff) force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); *damp=ob->pd->pdef_sbdamp; if (facedist > 0.0f) { *damp *= (1.0f - facedist/outerfacethickness); Vec3PlusStVec(outerforceaccu, force_mag_norm, d_nvect); deflected = 3; } else { Vec3PlusStVec(innerforceaccu, force_mag_norm, d_nvect); if (deflected < 2) deflected = 2; } if ((mprevvert) && (*damp > 0.0f)) { choose_winner(ve, opco, nv1, nv2, nv3, vv1, vv2, vv3); add_v3_v3(avel, ve); cavel ++; } *intrusion += facedist; ci++; } } if (mface->v4) { /* quad */ /* switch origin to be nv4 */ sub_v3_v3v3(edge1, nv3, nv4); sub_v3_v3v3(edge2, nv1, nv4); sub_v3_v3v3(dv1, opco, nv4); /* abuse dv1 to have vertex in question at *origin* of triangle */ cross_v3_v3v3(d_nvect, edge2, edge1); /* n_mag = */ /* UNUSED */ normalize_v3(d_nvect); facedist = dot_v3v3(dv1, d_nvect); if ((facedist > innerfacethickness) && (facedist < outerfacethickness)) { if (isect_point_tri_prism_v3(opco, nv1, nv3, nv4) ) { force_mag_norm =(float)exp(-ee*facedist); if (facedist > outerfacethickness*ff) force_mag_norm =(float)force_mag_norm*fa*(facedist - outerfacethickness)*(facedist - outerfacethickness); *damp=ob->pd->pdef_sbdamp; if (facedist > 0.0f) { *damp *= (1.0f - facedist/outerfacethickness); Vec3PlusStVec(outerforceaccu, force_mag_norm, d_nvect); deflected = 3; } else { Vec3PlusStVec(innerforceaccu, force_mag_norm, d_nvect); if (deflected < 2) deflected = 2; } if ((mprevvert) && (*damp > 0.0f)) { choose_winner(ve, opco, nv1, nv3, nv4, vv1, vv3, vv4); add_v3_v3(avel, ve); cavel ++; } *intrusion += facedist; ci++; } } if ((deflected < 2)&& (G.debug_value != 444)) { /* we did not hit a face until now */ /* see if 'outer' hits an edge */ float dist; closest_to_line_segment_v3(ve, opco, nv1, nv2); sub_v3_v3v3(ve, opco, ve); dist = normalize_v3(ve); if ((dist < outerfacethickness)&&(dist < mindistedge )) { copy_v3_v3(coledge, ve); mindistedge = dist, deflected=1; } closest_to_line_segment_v3(ve, opco, nv2, nv3); sub_v3_v3v3(ve, opco, ve); dist = normalize_v3(ve); if ((dist < outerfacethickness)&&(dist < mindistedge )) { copy_v3_v3(coledge, ve); mindistedge = dist, deflected=1; } closest_to_line_segment_v3(ve, opco, nv3, nv1); sub_v3_v3v3(ve, opco, ve); dist = normalize_v3(ve); if ((dist < outerfacethickness)&&(dist < mindistedge )) { copy_v3_v3(coledge, ve); mindistedge = dist, deflected=1; } if (mface->v4) { /* quad */ closest_to_line_segment_v3(ve, opco, nv3, nv4); sub_v3_v3v3(ve, opco, ve); dist = normalize_v3(ve); if ((dist < outerfacethickness)&&(dist < mindistedge )) { copy_v3_v3(coledge, ve); mindistedge = dist, deflected=1; } closest_to_line_segment_v3(ve, opco, nv1, nv4); sub_v3_v3v3(ve, opco, ve); dist = normalize_v3(ve); if ((dist < outerfacethickness)&&(dist < mindistedge )) { copy_v3_v3(coledge, ve); mindistedge = dist, deflected=1; } } } } mface++; mima++; }/* while a */ } /* if (ob->pd && ob->pd->deflect) */ BLI_ghashIterator_step(ihash); } } /* while () */ if (deflected == 1) { // no face but 'outer' edge cylinder sees vert force_mag_norm =(float)exp(-ee*mindistedge); if (mindistedge > outerfacethickness*ff) force_mag_norm =(float)force_mag_norm*fa*(mindistedge - outerfacethickness)*(mindistedge - outerfacethickness); Vec3PlusStVec(force, force_mag_norm, coledge); *damp=ob->pd->pdef_sbdamp; if (mindistedge > 0.0f) { *damp *= (1.0f - mindistedge/outerfacethickness); } } if (deflected == 2) { // face inner detected add_v3_v3(force, innerforceaccu); } if (deflected == 3) { // face outer detected add_v3_v3(force, outerforceaccu); } BLI_ghashIterator_free(ihash); if (cavel) mul_v3_fl(avel, 1.0f/(float)cavel); copy_v3_v3(vel, avel); if (ci) *intrusion /= ci; if (deflected) { normalize_v3_v3(facenormal, force); } return deflected; } /* sandbox to plug in various deflection algos */ static int sb_deflect_face(Object *ob, float *actpos, float *facenormal, float *force, float *cf, float time, float *vel, float *intrusion) { float s_actpos[3]; int deflected; copy_v3_v3(s_actpos, actpos); deflected= sb_detect_vertex_collisionCached(s_actpos, facenormal, cf, force, ob->lay, ob, time, vel, intrusion); //deflected= sb_detect_vertex_collisionCachedEx(s_actpos, facenormal, cf, force, ob->lay, ob, time, vel, intrusion); return(deflected); } /* hiding this for now .. but the jacobian may pop up on other tasks .. so i'd like to keep it static void dfdx_spring(int ia, int ic, int op, float dir[3], float L, float len, float factor) { float m, delta_ij; int i, j; if (L < len) { for (i=0;i<3;i++) for (j=0;j<3;j++) { delta_ij = (i==j ? (1.0f): (0.0f)); m=factor*(dir[i]*dir[j] + (1-L/len)*(delta_ij - dir[i]*dir[j])); nlMatrixAdd(ia+i, op+ic+j, m); } } else { for (i=0;i<3;i++) for (j=0;j<3;j++) { m=factor*dir[i]*dir[j]; nlMatrixAdd(ia+i, op+ic+j, m); } } } static void dfdx_goal(int ia, int ic, int op, float factor) { int i; for (i=0;i<3;i++) nlMatrixAdd(ia+i, op+ic+i, factor); } static void dfdv_goal(int ia, int ic, float factor) { int i; for (i=0;i<3;i++) nlMatrixAdd(ia+i, ic+i, factor); } */ static void sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime), int nl_flags) { SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp1, *bp2; float dir[3], dvel[3]; float distance, forcefactor, kd, absvel, projvel, kw; #if 0 /* UNUSED */ int ia, ic; #endif /* prepare depending on which side of the spring we are on */ if (bpi == bs->v1) { bp1 = &sb->bpoint[bs->v1]; bp2 = &sb->bpoint[bs->v2]; #if 0 /* UNUSED */ ia =3*bs->v1; ic =3*bs->v2; #endif } else if (bpi == bs->v2) { bp1 = &sb->bpoint[bs->v2]; bp2 = &sb->bpoint[bs->v1]; #if 0 /* UNUSED */ ia =3*bs->v2; ic =3*bs->v1; #endif } else { /* TODO make this debug option */ /**/ printf("bodypoint is not attached to spring <*bs> --> sb_spring_force()\n"); return; } /* do bp1 <--> bp2 elastic */ sub_v3_v3v3(dir, bp1->pos, bp2->pos); distance = normalize_v3(dir); if (bs->len < distance) iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ else iks = 1.0f/(1.0f-sb->inpush)-1.0f ;/* inner spring constants function */ if (bs->len > 0.0f) /* check for degenerated springs */ forcefactor = iks/bs->len; else forcefactor = iks; kw = (bp1->springweight+bp2->springweight)/2.0f; kw = kw * kw; kw = kw * kw; switch (bs->springtype) { case SB_EDGE: case SB_HANDLE: forcefactor *= kw; break; case SB_BEND: forcefactor *=sb->secondspring*kw; break; case SB_STIFFQUAD: forcefactor *=sb->shearstiff*sb->shearstiff* kw; break; default: break; } Vec3PlusStVec(bp1->force, (bs->len - distance)*forcefactor, dir); /* do bp1 <--> bp2 viscous */ sub_v3_v3v3(dvel, bp1->vec, bp2->vec); kd = sb->infrict * sb_fric_force_scale(ob); absvel = normalize_v3(dvel); projvel = dot_v3v3(dir, dvel); kd *= absvel * projvel; Vec3PlusStVec(bp1->force, -kd, dir); /* do jacobian stuff if needed */ if (nl_flags & NLF_BUILD) { //int op =3*sb->totpoint; //float mvel = -forcetime*kd; //float mpos = -forcetime*forcefactor; /* depending on my pos */ // dfdx_spring(ia, ia, op, dir, bs->len, distance, -mpos); /* depending on my vel */ // dfdv_goal(ia, ia, mvel); // well that ignores geometie if (bp2->goal < SOFTGOALSNAP) { /* omit this bp when it snaps */ /* depending on other pos */ // dfdx_spring(ia, ic, op, dir, bs->len, distance, mpos); /* depending on other vel */ // dfdv_goal(ia, ia, -mvel); // well that ignores geometie } } } /* since this is definitely the most CPU consuming task here .. try to spread it */ /* core function _softbody_calc_forces_slice_in_a_thread */ /* result is int to be able to flag user break */ 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) { float iks; int bb, do_selfcollision, do_springcollision, do_aero; int number_of_points_here = ilast - ifirst; SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; /* intitialize */ if (sb) { /* check conditions for various options */ /* +++ could be done on object level to squeeze out the last bits of it */ do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); /* --- could be done on object level to squeeze out the last bits of it */ } else { printf("Error expected a SB here\n"); return (999); } /* debugerin */ if (sb->totpoint < ifirst) { printf("Aye 998"); return (998); } /* debugerin */ bp = &sb->bpoint[ifirst]; for (bb=number_of_points_here; bb>0; bb--, bp++) { /* clear forces accumulator */ bp->force[0] = bp->force[1] = bp->force[2] = 0.0; /* naive ball self collision */ /* needs to be done if goal snaps or not */ if (do_selfcollision) { int attached; BodyPoint *obp; BodySpring *bs; int c, b; float velcenter[3], dvel[3], def[3]; float distance; float compare; float bstune = sb->ballstiff; for (c=sb->totpoint, obp= sb->bpoint; c>=ifirst+bb; c--, obp++) { compare = (obp->colball + bp->colball); sub_v3_v3v3(def, bp->pos, obp->pos); /* rather check the AABBoxes before ever calulating the real distance */ /* mathematically it is completly nuts, but performance is pretty much (3) times faster */ if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; distance = normalize_v3(def); if (distance < compare ) { /* exclude body points attached with a spring */ attached = 0; for (b=obp->nofsprings;b>0;b--) { bs = sb->bspring + obp->springs[b-1]; if (( ilast-bb == bs->v2) || ( ilast-bb == bs->v1)) { attached=1; continue;} } if (!attached) { float f = bstune / (distance) + bstune / (compare * compare) * distance - 2.0f * bstune / compare; mid_v3_v3v3(velcenter, bp->vec, obp->vec); sub_v3_v3v3(dvel, velcenter, bp->vec); mul_v3_fl(dvel, _final_mass(ob, bp)); Vec3PlusStVec(bp->force, f*(1.0f-sb->balldamp), def); Vec3PlusStVec(bp->force, sb->balldamp, dvel); /* exploit force(a, b) == -force(b, a) part2/2 */ sub_v3_v3v3(dvel, velcenter, obp->vec); mul_v3_fl(dvel, _final_mass(ob, bp)); Vec3PlusStVec(obp->force, sb->balldamp, dvel); Vec3PlusStVec(obp->force, -f*(1.0f-sb->balldamp), def); } } } } /* naive ball self collision done */ if (_final_goal(ob, bp) < SOFTGOALSNAP) { /* omit this bp when it snaps */ float auxvect[3]; float velgoal[3]; /* do goal stuff */ if (ob->softflag & OB_SB_GOAL) { /* true elastic goal */ float ks, kd; sub_v3_v3v3(auxvect, bp->pos, bp->origT); ks = 1.0f / (1.0f - _final_goal(ob, bp) * sb->goalspring) - 1.0f; bp->force[0]+= -ks*(auxvect[0]); bp->force[1]+= -ks*(auxvect[1]); bp->force[2]+= -ks*(auxvect[2]); /* calulate damping forces generated by goals*/ sub_v3_v3v3(velgoal, bp->origS, bp->origE); kd = sb->goalfrict * sb_fric_force_scale(ob); add_v3_v3v3(auxvect, velgoal, bp->vec); if (forcetime > 0.0f) { /* make sure friction does not become rocket motor on time reversal */ bp->force[0]-= kd * (auxvect[0]); bp->force[1]-= kd * (auxvect[1]); bp->force[2]-= kd * (auxvect[2]); } else { bp->force[0]-= kd * (velgoal[0] - bp->vec[0]); bp->force[1]-= kd * (velgoal[1] - bp->vec[1]); bp->force[2]-= kd * (velgoal[2] - bp->vec[2]); } } /* done goal stuff */ /* gravitation */ if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) { float gravity[3]; copy_v3_v3(gravity, scene->physics_settings.gravity); mul_v3_fl(gravity, sb_grav_force_scale(ob)*_final_mass(ob, bp)*sb->effector_weights->global_gravity); /* individual mass of node here */ add_v3_v3(bp->force, gravity); } /* particle field & vortex */ if (do_effector) { EffectedPoint epoint; float kd; float force[3] = {0.0f, 0.0f, 0.0f}; float speed[3] = {0.0f, 0.0f, 0.0f}; float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint); pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); /* apply forcefield*/ mul_v3_fl(force, fieldfactor* eval_sb_fric_force_scale); add_v3_v3(bp->force, force); /* BP friction in moving media */ kd= sb->mediafrict* eval_sb_fric_force_scale; bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale); bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale); bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale); /* now we'll have nice centrifugal effect for vortex */ } else { /* BP friction in media (not) moving*/ float kd = sb->mediafrict* sb_fric_force_scale(ob); /* assume it to be proportional to actual velocity */ bp->force[0]-= bp->vec[0]*kd; bp->force[1]-= bp->vec[1]*kd; bp->force[2]-= bp->vec[2]*kd; /* friction in media done */ } /* +++cached collision targets */ bp->choke = 0.0f; bp->choke2 = 0.0f; bp->loc_flag &= ~SBF_DOFUZZY; if (do_deflector && !(bp->loc_flag & SBF_OUTOFCOLLISION) ) { 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; float kd = 1.0f; if (sb_deflect_face(ob, bp->pos, facenormal, defforce, &cf, timenow, vel, &intrusion)) { if (intrusion < 0.0f) { sb->scratch->flag |= SBF_DOFUZZY; bp->loc_flag |= SBF_DOFUZZY; bp->choke = sb->choke*0.01f; } sub_v3_v3v3(cfforce, bp->vec, vel); Vec3PlusStVec(bp->force, -cf*50.0f, cfforce); Vec3PlusStVec(bp->force, kd, defforce); } } /* ---cached collision targets */ /* +++springs */ iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ if (ob->softflag & OB_SB_EDGES) { if (sb->bspring) { /* spring list exists at all ? */ int b; BodySpring *bs; for (b=bp->nofsprings;b>0;b--) { bs = sb->bspring + bp->springs[b-1]; if (do_springcollision || do_aero) { add_v3_v3(bp->force, bs->ext_force); if (bs->flag & BSF_INTERSECT) bp->choke = bs->cf; } // sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float forcetime, int nl_flags) sb_spring_force(ob, ilast-bb, bs, iks, forcetime, 0); }/* loop springs */ }/* existing spring list */ }/*any edges*/ /* ---springs */ }/*omit on snap */ }/*loop all bp's*/ return 0; /*done fine*/ } static void *exec_softbody_calc_forces(void *data) { SB_thread_context *pctx = (SB_thread_context*)data; _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); return NULL; } 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) { ListBase threads; SB_thread_context *sb_threads; int i, totthread, left, dec; int lowpoints =100; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */ /* figure the number of threads while preventing pretty pointless threading overhead */ totthread= BKE_scene_num_threads(scene); /* what if we got zillions of CPUs running but less to spread*/ while ((totpoint/totthread < lowpoints) && (totthread > 1)) { totthread--; } /* printf("sb_cf_threads_run spawning %d threads\n", totthread); */ sb_threads= MEM_callocN(sizeof(SB_thread_context)*totthread, "SBThread"); memset(sb_threads, 0, sizeof(SB_thread_context)*totthread); left = totpoint; dec = totpoint/totthread +1; for (i=0; i0) { sb_threads[i].ifirst = left; } else sb_threads[i].ifirst = 0; sb_threads[i].do_effector = do_effector; sb_threads[i].do_deflector = do_deflector; sb_threads[i].fieldfactor = fieldfactor; sb_threads[i].windfactor = windfactor; sb_threads[i].nr= i; sb_threads[i].tot= totthread; } if (totthread > 1) { BLI_init_threads(&threads, exec_softbody_calc_forces, totthread); for (i=0; ivec bp->pos in here ! * this will ruin adaptive stepsize AKA heun! (BM) */ SoftBody *sb= ob->soft; /* is supposed to be there */ /*BodyPoint *bproot;*/ /* UNUSED */ ListBase *do_effector = NULL; /* float gravity; */ /* UNUSED */ /* float iks; */ float fieldfactor = -1.0f, windfactor = 0.25; int do_deflector /*, do_selfcollision*/, do_springcollision, do_aero; /* gravity = sb->grav * sb_grav_force_scale(ob); */ /* UNUSED */ /* check conditions for various options */ do_deflector= query_external_colliders(scene, ob); /* do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); */ /* UNUSED */ do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); /* iks = 1.0f/(1.0f-sb->inspring)-1.0f; */ /* inner spring constants function */ /* UNUSED */ /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */ if (do_springcollision || do_aero) sb_sfesf_threads_run(scene, ob, timenow, sb->totspring, NULL); /* after spring scan because it uses Effoctors too */ do_effector= pdInitEffectors(scene, ob, NULL, sb->effector_weights, true); if (do_deflector) { float defforce[3]; do_deflector = sb_detect_aabb_collisionCached(defforce, ob->lay, ob, timenow); } sb_cf_threads_run(scene, ob, forcetime, timenow, sb->totpoint, NULL, do_effector, do_deflector, fieldfactor, windfactor); /* finally add forces caused by face collision */ if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob, timenow); /* finish matrix and solve */ pdEndEffectors(&do_effector); } static void softbody_calc_forces(Scene *scene, Object *ob, float forcetime, float timenow, int nl_flags) { /* redirection to the new threaded Version */ if (!(G.debug_value & 0x10)) { // 16 softbody_calc_forcesEx(scene, ob, forcetime, timenow, nl_flags); return; } else { /* so the following will die */ /* |||||||||||||||||||||||||| */ /* VVVVVVVVVVVVVVVVVVVVVVVVVV */ /*backward compatibility note: fixing bug [17428] which forces adaptive step size to tiny steps in some situations .. keeping G.debug_value==17 0x11 option for old files 'needing' the bug*/ /* rule we never alter free variables :bp->vec bp->pos in here ! * this will ruin adaptive stepsize AKA heun! (BM) */ SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; /* BodyPoint *bproot; */ /* UNUSED */ BodySpring *bs; ListBase *do_effector = NULL; float iks, ks, kd, gravity[3] = {0.0f, 0.0f, 0.0f}; float fieldfactor = -1.0f, windfactor = 0.25f; float tune = sb->ballstiff; int a, b, do_deflector, do_selfcollision, do_springcollision, do_aero; /* jacobian NLboolean success; if (nl_flags) { nlBegin(NL_SYSTEM); nlBegin(NL_MATRIX); } */ if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) { copy_v3_v3(gravity, scene->physics_settings.gravity); mul_v3_fl(gravity, sb_grav_force_scale(ob)*sb->effector_weights->global_gravity); } /* check conditions for various options */ do_deflector= query_external_colliders(scene, ob); do_selfcollision=((ob->softflag & OB_SB_EDGES) && (sb->bspring)&& (ob->softflag & OB_SB_SELF)); do_springcollision=do_deflector && (ob->softflag & OB_SB_EDGES) &&(ob->softflag & OB_SB_EDGECOLL); do_aero=((sb->aeroedge)&& (ob->softflag & OB_SB_EDGES)); iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */ /* bproot= sb->bpoint; */ /* need this for proper spring addressing */ /* UNUSED */ if (do_springcollision || do_aero) scan_for_ext_spring_forces(scene, ob, timenow); /* after spring scan because it uses Effoctors too */ do_effector= pdInitEffectors(scene, ob, NULL, ob->soft->effector_weights, true); if (do_deflector) { float defforce[3]; do_deflector = sb_detect_aabb_collisionCached(defforce, ob->lay, ob, timenow); } for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { /* clear forces accumulator */ bp->force[0] = bp->force[1] = bp->force[2] = 0.0; if (nl_flags & NLF_BUILD) { //int ia =3*(sb->totpoint-a); //int op =3*sb->totpoint; /* dF/dV = v */ /* jacobioan nlMatrixAdd(op+ia, ia, -forcetime); nlMatrixAdd(op+ia+1, ia+1, -forcetime); nlMatrixAdd(op+ia+2, ia+2, -forcetime); nlMatrixAdd(ia, ia, 1); nlMatrixAdd(ia+1, ia+1, 1); nlMatrixAdd(ia+2, ia+2, 1); nlMatrixAdd(op+ia, op+ia, 1); nlMatrixAdd(op+ia+1, op+ia+1, 1); nlMatrixAdd(op+ia+2, op+ia+2, 1); */ } /* naive ball self collision */ /* needs to be done if goal snaps or not */ if (do_selfcollision) { int attached; BodyPoint *obp; int c, b; float velcenter[3], dvel[3], def[3]; float distance; float compare; for (c=sb->totpoint, obp= sb->bpoint; c>=a; c--, obp++) { //if ((bp->octantflag & obp->octantflag) == 0) continue; compare = (obp->colball + bp->colball); sub_v3_v3v3(def, bp->pos, obp->pos); /* rather check the AABBoxes before ever calulating the real distance */ /* mathematically it is completly nuts, but performance is pretty much (3) times faster */ if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue; distance = normalize_v3(def); if (distance < compare ) { /* exclude body points attached with a spring */ attached = 0; for (b=obp->nofsprings;b>0;b--) { bs = sb->bspring + obp->springs[b-1]; if (( sb->totpoint-a == bs->v2) || ( sb->totpoint-a == bs->v1)) { attached=1; continue;} } if (!attached) { float f = tune / (distance) + tune / (compare * compare) * distance - 2.0f * tune/compare; mid_v3_v3v3(velcenter, bp->vec, obp->vec); sub_v3_v3v3(dvel, velcenter, bp->vec); mul_v3_fl(dvel, _final_mass(ob, bp)); Vec3PlusStVec(bp->force, f*(1.0f-sb->balldamp), def); Vec3PlusStVec(bp->force, sb->balldamp, dvel); if (nl_flags & NLF_BUILD) { //int ia =3*(sb->totpoint-a); //int ic =3*(sb->totpoint-c); //int op =3*sb->totpoint; //float mvel = forcetime*sb->nodemass*sb->balldamp; //float mpos = forcetime*tune*(1.0f-sb->balldamp); /*some quick and dirty entries to the jacobian*/ //dfdx_goal(ia, ia, op, mpos); //dfdv_goal(ia, ia, mvel); /* exploit force(a, b) == -force(b, a) part1/2 */ //dfdx_goal(ic, ic, op, mpos); //dfdv_goal(ic, ic, mvel); /*TODO sit down an X-out the true jacobian entries*/ /*well does not make to much sense because the eigenvalues of the jacobian go negative; and negative eigenvalues on a complex iterative system z(n+1)=A * z(n) give imaginary roots in the charcateristic polynom --> solutions that to z(t)=u(t)* exp ( i omega t) --> oscilations we don't want here where u(t) is a unknown amplitude function (worst case rising fast) */ } /* exploit force(a, b) == -force(b, a) part2/2 */ sub_v3_v3v3(dvel, velcenter, obp->vec); mul_v3_fl(dvel, (_final_mass(ob, bp)+_final_mass(ob, obp))/2.0f); Vec3PlusStVec(obp->force, sb->balldamp, dvel); Vec3PlusStVec(obp->force, -f*(1.0f-sb->balldamp), def); } } } } /* naive ball self collision done */ if (_final_goal(ob, bp) < SOFTGOALSNAP) { /* omit this bp when it snaps */ float auxvect[3]; float velgoal[3]; /* do goal stuff */ if (ob->softflag & OB_SB_GOAL) { /* true elastic goal */ sub_v3_v3v3(auxvect, bp->pos, bp->origT); ks = 1.0f / (1.0f- _final_goal(ob, bp) * sb->goalspring) - 1.0f; bp->force[0]+= -ks*(auxvect[0]); bp->force[1]+= -ks*(auxvect[1]); bp->force[2]+= -ks*(auxvect[2]); if (nl_flags & NLF_BUILD) { //int ia =3*(sb->totpoint-a); //int op =3*(sb->totpoint); /* depending on my pos */ //dfdx_goal(ia, ia, op, ks*forcetime); } /* calulate damping forces generated by goals*/ sub_v3_v3v3(velgoal, bp->origS, bp->origE); kd = sb->goalfrict * sb_fric_force_scale(ob); add_v3_v3v3(auxvect, velgoal, bp->vec); if (forcetime > 0.0f) { /* make sure friction does not become rocket motor on time reversal */ bp->force[0]-= kd * (auxvect[0]); bp->force[1]-= kd * (auxvect[1]); bp->force[2]-= kd * (auxvect[2]); if (nl_flags & NLF_BUILD) { //int ia =3*(sb->totpoint-a); normalize_v3(auxvect); /* depending on my vel */ //dfdv_goal(ia, ia, kd*forcetime); } } else { bp->force[0]-= kd * (velgoal[0] - bp->vec[0]); bp->force[1]-= kd * (velgoal[1] - bp->vec[1]); bp->force[2]-= kd * (velgoal[2] - bp->vec[2]); } } /* done goal stuff */ /* gravitation */ madd_v3_v3fl(bp->force, gravity, _final_mass(ob, bp)); /* individual mass of node here */ /* particle field & vortex */ if (do_effector) { EffectedPoint epoint; float force[3] = {0.0f, 0.0f, 0.0f}; float speed[3] = {0.0f, 0.0f, 0.0f}; float eval_sb_fric_force_scale = sb_fric_force_scale(ob); /* just for calling function once */ pd_point_from_soft(scene, bp->pos, bp->vec, sb->bpoint-bp, &epoint); pdDoEffectors(do_effector, NULL, sb->effector_weights, &epoint, force, speed); /* apply forcefield*/ mul_v3_fl(force, fieldfactor* eval_sb_fric_force_scale); add_v3_v3(bp->force, force); /* BP friction in moving media */ kd= sb->mediafrict* eval_sb_fric_force_scale; bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale); bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale); bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale); /* now we'll have nice centrifugal effect for vortex */ } else { /* BP friction in media (not) moving*/ kd= sb->mediafrict* sb_fric_force_scale(ob); /* assume it to be proportional to actual velocity */ bp->force[0]-= bp->vec[0]*kd; bp->force[1]-= bp->vec[1]*kd; bp->force[2]-= bp->vec[2]*kd; /* friction in media done */ if (nl_flags & NLF_BUILD) { //int ia =3*(sb->totpoint-a); /* da/dv = */ // nlMatrixAdd(ia, ia, forcetime*kd); // nlMatrixAdd(ia+1, ia+1, forcetime*kd); // nlMatrixAdd(ia+2, ia+2, forcetime*kd); } } /* +++cached collision targets */ bp->choke = 0.0f; bp->choke2 = 0.0f; bp->loc_flag &= ~SBF_DOFUZZY; if (do_deflector) { 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; kd = 1.0f; if (sb_deflect_face(ob, bp->pos, facenormal, defforce, &cf, timenow, vel, &intrusion)) { if ((!nl_flags)&&(intrusion < 0.0f)) { if (G.debug_value & 0x01) { // 17 we did check for bit 0x10 before /*fixing bug [17428] this forces adaptive step size to tiny steps in some situations .. keeping G.debug_value==17 option for old files 'needing' the bug */ /*bjornmose: uugh.. what an evil hack violation of the 'don't touch bp->pos in here' rule but works nice, like this--> we predict the solution being out of the collider in heun step No1 and leave the heun step No2 adapt to it so we kind of introduced a implicit solver for this case */ Vec3PlusStVec(bp->pos, -intrusion, facenormal); } else { sub_v3_v3v3(cfforce, bp->vec, vel); Vec3PlusStVec(bp->force, -cf*50.0f, cfforce); } sb->scratch->flag |= SBF_DOFUZZY; bp->loc_flag |= SBF_DOFUZZY; bp->choke = sb->choke*0.01f; } else { sub_v3_v3v3(cfforce, bp->vec, vel); Vec3PlusStVec(bp->force, -cf*50.0f, cfforce); } Vec3PlusStVec(bp->force, kd, defforce); if (nl_flags & NLF_BUILD) { // int ia =3*(sb->totpoint-a); // int op =3*sb->totpoint; //dfdx_goal(ia, ia, op, mpos); // don't do unless you know //dfdv_goal(ia, ia, -cf); } } } /* ---cached collision targets */ /* +++springs */ if (ob->softflag & OB_SB_EDGES) { if (sb->bspring) { /* spring list exists at all ? */ for (b=bp->nofsprings;b>0;b--) { bs = sb->bspring + bp->springs[b-1]; if (do_springcollision || do_aero) { add_v3_v3(bp->force, bs->ext_force); if (bs->flag & BSF_INTERSECT) bp->choke = bs->cf; } // sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float forcetime, int nl_flags) // rather remove nl_falgs from code .. will make things a lot cleaner sb_spring_force(ob, sb->totpoint-a, bs, iks, forcetime, 0); }/* loop springs */ }/* existing spring list */ }/*any edges*/ /* ---springs */ }/*omit on snap */ }/*loop all bp's*/ /* finally add forces caused by face collision */ if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob, timenow); /* finish matrix and solve */ #if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be useful .. so keep that BM if (nl_flags & NLF_SOLVE) { //double sct, sst=PIL_check_seconds_timer(); for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { int iv =3*(sb->totpoint-a); int ip =3*(2*sb->totpoint-a); int n; for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);} for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);} } nlEnd(NL_MATRIX); nlEnd(NL_SYSTEM); if ((G.debug_value == 32) && (nl_flags & NLF_BUILD)) { printf("####MEE#####\n"); nlPrintMatrix(); } success= nlSolveAdvanced(NULL, 1); // nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */ if (success) { float f; int index =0; /* for debug purpose .. anyhow cropping B vector looks like working */ if (G.debug_value ==32) for (a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { f=nlGetVariable(0, index); printf("(%f ", f);index++; f=nlGetVariable(0, index); printf("%f ", f);index++; f=nlGetVariable(0, index); printf("%f)", f);index++; } index =0; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { f=nlGetVariable(0, index); bp->impdv[0] = f; index++; f=nlGetVariable(0, index); bp->impdv[1] = f; index++; f=nlGetVariable(0, index); bp->impdv[2] = f; index++; } /* for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { f=nlGetVariable(0, index); bp->impdx[0] = f; index++; f=nlGetVariable(0, index); bp->impdx[1] = f; index++; f=nlGetVariable(0, index); bp->impdx[2] = f; index++; } */ } else { printf("Matrix inversion failed\n"); for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(bp->impdv, bp->force); } } //sct=PIL_check_seconds_timer(); //if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r", sct-sst, ob->id.name); } /* cleanup */ #endif pdEndEffectors(&do_effector); } } static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags) { /* time evolution */ /* actually does an explicit euler step mode == 0 */ /* or heun ~ 2nd order runge-kutta steps, mode 1, 2 */ SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; float dx[3] = {0}, dv[3], aabbmin[3], aabbmax[3], cm[3] = {0.0f, 0.0f, 0.0f}; float timeovermass/*, freezeloc=0.00001f, freezeforce=0.00000000001f*/; float maxerrpos= 0.0f, maxerrvel = 0.0f; int a, fuzzy=0; forcetime *= sb_time_scale(ob); aabbmin[0]=aabbmin[1]=aabbmin[2] = 1e20f; aabbmax[0]=aabbmax[1]=aabbmax[2] = -1e20f; /* old one with homogeneous masses */ /* claim a minimum mass for vertex */ /* if (sb->nodemass > 0.009999f) timeovermass = forcetime/sb->nodemass; else timeovermass = forcetime/0.009999f; */ for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { /* now we have individual masses */ /* claim a minimum mass for vertex */ if (_final_mass(ob, bp) > 0.009999f) timeovermass = forcetime/_final_mass(ob, bp); else timeovermass = forcetime/0.009999f; if (_final_goal(ob, bp) < SOFTGOALSNAP) { /* this makes t~ = t */ if (mid_flags & MID_PRESERVE) copy_v3_v3(dx, bp->vec); /* so here is (v)' = a(cceleration) = sum(F_springs)/m + gravitation + some friction forces + more forces*/ /* the ( ... )' operator denotes derivate respective time */ /* the euler step for velocity then becomes */ /* v(t + dt) = v(t) + a(t) * dt */ mul_v3_fl(bp->force, timeovermass);/* individual mass of node here */ /* some nasty if's to have heun in here too */ copy_v3_v3(dv, bp->force); if (mode == 1) { copy_v3_v3(bp->prevvec, bp->vec); copy_v3_v3(bp->prevdv, dv); } if (mode ==2) { /* be optimistic and execute step */ bp->vec[0] = bp->prevvec[0] + 0.5f * (dv[0] + bp->prevdv[0]); bp->vec[1] = bp->prevvec[1] + 0.5f * (dv[1] + bp->prevdv[1]); bp->vec[2] = bp->prevvec[2] + 0.5f * (dv[2] + bp->prevdv[2]); /* compare euler to heun to estimate error for step sizing */ maxerrvel = max_ff(maxerrvel, fabsf(dv[0] - bp->prevdv[0])); maxerrvel = max_ff(maxerrvel, fabsf(dv[1] - bp->prevdv[1])); maxerrvel = max_ff(maxerrvel, fabsf(dv[2] - bp->prevdv[2])); } else { add_v3_v3(bp->vec, bp->force); } /* this makes t~ = t+dt */ if (!(mid_flags & MID_PRESERVE)) copy_v3_v3(dx, bp->vec); /* so here is (x)'= v(elocity) */ /* the euler step for location then becomes */ /* x(t + dt) = x(t) + v(t~) * dt */ mul_v3_fl(dx, forcetime); /* the freezer coming sooner or later */ #if 0 if ((dot_v3v3(dx, dx)force, bp->force)frozen /=2; } else { bp->frozen = min_ff(bp->frozen*1.05f, 1.0f); } mul_v3_fl(dx, bp->frozen); #endif /* again some nasty if's to have heun in here too */ if (mode ==1) { copy_v3_v3(bp->prevpos, bp->pos); copy_v3_v3(bp->prevdx, dx); } if (mode ==2) { bp->pos[0] = bp->prevpos[0] + 0.5f * ( dx[0] + bp->prevdx[0]); bp->pos[1] = bp->prevpos[1] + 0.5f * ( dx[1] + bp->prevdx[1]); bp->pos[2] = bp->prevpos[2] + 0.5f * ( dx[2] + bp->prevdx[2]); maxerrpos = max_ff(maxerrpos, fabsf(dx[0] - bp->prevdx[0])); maxerrpos = max_ff(maxerrpos, fabsf(dx[1] - bp->prevdx[1])); maxerrpos = max_ff(maxerrpos, fabsf(dx[2] - bp->prevdx[2])); /* bp->choke is set when we need to pull a vertex or edge out of the collider. * the collider object signals to get out by pushing hard. on the other hand * we don't want to end up in deep space so we add some * to balance that out */ if (bp->choke2 > 0.0f) { mul_v3_fl(bp->vec, (1.0f - bp->choke2)); } if (bp->choke > 0.0f) { mul_v3_fl(bp->vec, (1.0f - bp->choke)); } } else { add_v3_v3(bp->pos, dx);} }/*snap*/ /* so while we are looping BPs anyway do statistics on the fly */ minmax_v3v3_v3(aabbmin, aabbmax, bp->pos); if (bp->loc_flag & SBF_DOFUZZY) fuzzy =1; } /*for*/ if (sb->totpoint) mul_v3_fl(cm, 1.0f/sb->totpoint); if (sb->scratch) { copy_v3_v3(sb->scratch->aabbmin, aabbmin); copy_v3_v3(sb->scratch->aabbmax, aabbmax); } if (err) { /* so step size will be controlled by biggest difference in slope */ if (sb->solverflags & SBSO_OLDERR) *err = max_ff(maxerrpos, maxerrvel); else *err = maxerrpos; //printf("EP %f EV %f\n", maxerrpos, maxerrvel); if (fuzzy) { *err /= sb->fuzzyness; } } } /* used by heun when it overshoots */ static void softbody_restore_prev_step(Object *ob) { SoftBody *sb= ob->soft; /* is supposed to be there*/ BodyPoint *bp; int a; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(bp->vec, bp->prevvec); copy_v3_v3(bp->pos, bp->prevpos); } } #if 0 static void softbody_store_step(Object *ob) { SoftBody *sb= ob->soft; /* is supposed to be there*/ BodyPoint *bp; int a; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(bp->prevvec, bp->vec); copy_v3_v3(bp->prevpos, bp->pos); } } /* used by predictors and correctors */ static void softbody_store_state(Object *ob, float *ppos, float *pvel) { SoftBody *sb= ob->soft; /* is supposed to be there*/ BodyPoint *bp; int a; float *pp=ppos, *pv=pvel; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(pv, bp->vec); pv+=3; copy_v3_v3(pp, bp->pos); pp+=3; } } /* used by predictors and correctors */ static void softbody_retrieve_state(Object *ob, float *ppos, float *pvel) { SoftBody *sb= ob->soft; /* is supposed to be there*/ BodyPoint *bp; int a; float *pp=ppos, *pv=pvel; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(bp->vec, pv); pv+=3; copy_v3_v3(bp->pos, pp); pp+=3; } } /* used by predictors and correctors */ static void softbody_swap_state(Object *ob, float *ppos, float *pvel) { SoftBody *sb= ob->soft; /* is supposed to be there*/ BodyPoint *bp; int a; float *pp=ppos, *pv=pvel; float temp[3]; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { copy_v3_v3(temp, bp->vec); copy_v3_v3(bp->vec, pv); copy_v3_v3(pv, temp); pv+=3; copy_v3_v3(temp, bp->pos); copy_v3_v3(bp->pos, pp); copy_v3_v3(pp, temp); pp+=3; } } #endif /* care for bodypoints taken out of the 'ordinary' solver step ** because they are screwed to goal by bolts ** they just need to move along with the goal in time ** we need to adjust them on sub frame timing in solver ** so now when frame is done .. put 'em to the position at the end of frame */ static void softbody_apply_goalsnap(Object *ob) { SoftBody *sb= ob->soft; /* is supposed to be there */ BodyPoint *bp; int a; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { if (_final_goal(ob, bp) >= SOFTGOALSNAP) { copy_v3_v3(bp->prevpos, bp->pos); copy_v3_v3(bp->pos, bp->origT); } } } static void apply_spring_memory(Object *ob) { SoftBody *sb = ob->soft; BodySpring *bs; BodyPoint *bp1, *bp2; int a; float b, l, r; if (sb && sb->totspring) { b = sb->plastic; for (a=0; atotspring; a++) { bs = &sb->bspring[a]; bp1 =&sb->bpoint[bs->v1]; bp2 =&sb->bpoint[bs->v2]; l = len_v3v3(bp1->pos, bp2->pos); r = bs->len/l; if (( r > 1.05f) || (r < 0.95f)) { bs->len = ((100.0f - b) * bs->len + b*l)/100.0f; } } } } /* expects full initialized softbody */ static void interpolate_exciter(Object *ob, int timescale, int time) { SoftBody *sb= ob->soft; BodyPoint *bp; float f; int a; f = (float)time/(float)timescale; for (a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) { bp->origT[0] = bp->origS[0] + f*(bp->origE[0] - bp->origS[0]); bp->origT[1] = bp->origS[1] + f*(bp->origE[1] - bp->origS[1]); bp->origT[2] = bp->origS[2] + f*(bp->origE[2] - bp->origS[2]); if (_final_goal(ob, bp) >= SOFTGOALSNAP) { bp->vec[0] = bp->origE[0] - bp->origS[0]; bp->vec[1] = bp->origE[1] - bp->origS[1]; bp->vec[2] = bp->origE[2] - bp->origS[2]; } } } /* ************ convertors ********** */ /* for each object type we need; - xxxx_to_softbody(Object *ob) : a full (new) copy, creates SB geometry */ static void get_scalar_from_vertexgroup(Object *ob, int vertID, int groupindex, float *target) /* result 0 on success, else indicates error number -- kind of *inverse* result defintion, -- but this way we can signal error condition to caller -- and yes this function must not be here but in a *vertex group module* */ { MDeformVert *dv= NULL; int i; /* spot the vert in deform vert list at mesh */ if (ob->type==OB_MESH) { Mesh *me= ob->data; if (me->dvert) dv = me->dvert + vertID; } else if (ob->type==OB_LATTICE) { /* not yet supported in softbody btw */ Lattice *lt= ob->data; if (lt->dvert) dv = lt->dvert + vertID; } if (dv) { /* Lets see if this vert is in the weight group */ for (i=0; itotweight; i++) { if (dv->dw[i].def_nr == groupindex) { *target= dv->dw[i].weight; /* got it ! */ break; } } } } /* Resetting a Mesh SB object's springs */ /* Spring length are caculted from'raw' mesh vertices that are NOT altered by modifier stack. */ static void springs_from_mesh(Object *ob) { SoftBody *sb; Mesh *me= ob->data; BodyPoint *bp; int a; float scale =1.0f; sb= ob->soft; if (me && sb) { /* using bp->origS as a container for spring calcualtions here * will be overwritten sbObjectStep() to receive * actual modifier stack positions */ if (me->totvert) { bp= ob->soft->bpoint; for (a=0; atotvert; a++, bp++) { copy_v3_v3(bp->origS, me->mvert[a].co); mul_m4_v3(ob->obmat, bp->origS); } } /* recalculate spring length for meshes here */ /* public version shrink to fit */ if (sb->springpreload != 0 ) { scale = sb->springpreload / 100.0f; } for (a=0; atotspring; a++) { BodySpring *bs = &sb->bspring[a]; bs->len= scale*len_v3v3(sb->bpoint[bs->v1].origS, sb->bpoint[bs->v2].origS); } } } /* makes totally fresh start situation */ static void mesh_to_softbody(Scene *scene, Object *ob) { SoftBody *sb; Mesh *me= ob->data; MEdge *medge= me->medge; BodyPoint *bp; BodySpring *bs; int a, totedge; BKE_mesh_tessface_ensure(me); if (ob->softflag & OB_SB_EDGES) totedge= me->totedge; else totedge= 0; /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ renew_softbody(scene, ob, me->totvert, totedge); /* we always make body points */ sb= ob->soft; bp= sb->bpoint; for (a=0; atotvert; a++, bp++) { /* get scalar values needed *per vertex* from vertex group functions, * so we can *paint* them nicly .. * they are normalized [0.0..1.0] so may be we need amplitude for scale * which can be done by caller but still .. i'd like it to go this way */ if ((ob->softflag & OB_SB_GOAL) && sb->vertgroup) { /* even this is a deprecated evil hack */ /* I'd like to have it .. if (sb->namedVG_Goal[0]) */ get_scalar_from_vertexgroup(ob, a, sb->vertgroup - 1, &bp->goal); /* do this always, regardless successful read from vertex group */ /* this is where '2.5 every thing is animatable' goes wrong in the first place jow_go_for2_5 */ /* 1st coding action to take : move this to frame level */ /* reads: leave the bp->goal as it was read from vertex group / or default .. we will need it at per frame call */ /* should be fixed for meshes */ // bp->goal= sb->mingoal + bp->goal*goalfac; /* do not do here jow_go_for2_5 */ } else { /* in consequence if no group was set .. but we want to animate it laters */ /* logically attach to goal with default first */ if (ob->softflag & OB_SB_GOAL) {bp->goal = sb->defgoal;} } /* to proof the concept * this enables per vertex *mass painting* */ if (sb->namedVG_Mass[0]) { int defgrp_index = defgroup_name_index(ob, sb->namedVG_Mass); /* printf("VGN %s %d\n", sb->namedVG_Mass, defgrp_index); */ if (defgrp_index != -1) { get_scalar_from_vertexgroup(ob, a, defgrp_index, &bp->mass); /* 2.5 bp->mass = bp->mass * sb->nodemass; */ /* printf("bp->mass %f\n", bp->mass); */ } } /* first set the default */ bp->springweight = 1.0f; if (sb->namedVG_Spring_K[0]) { int defgrp_index = defgroup_name_index(ob, sb->namedVG_Spring_K); //printf("VGN %s %d\n", sb->namedVG_Spring_K, defgrp_index); if (defgrp_index != -1) { get_scalar_from_vertexgroup(ob, a, defgrp_index , &bp->springweight); //printf("bp->springweight %f\n", bp->springweight); } } } /* but we only optionally add body edge springs */ if (ob->softflag & OB_SB_EDGES) { if (medge) { bs= sb->bspring; for (a=me->totedge; a>0; a--, medge++, bs++) { bs->v1= medge->v1; bs->v2= medge->v2; bs->springtype=SB_EDGE; } /* insert *diagonal* springs in quads if desired */ if (ob->softflag & OB_SB_QUADS) { add_mesh_quad_diag_springs(ob); } build_bps_springlist(ob); /* scan for springs attached to bodypoints ONCE */ /* insert *other second order* springs if desired */ if (sb->secondspring > 0.0000001f) { add_2nd_order_springs(ob, sb->secondspring); /* exploits the first run of build_bps_springlist(ob);*/ build_bps_springlist(ob); /* yes we need to do it again*/ } springs_from_mesh(ob); /* write the 'rest'-length of the springs */ if (ob->softflag & OB_SB_SELF) {calculate_collision_balls(ob);} } } } static void mesh_faces_to_scratch(Object *ob) { SoftBody *sb= ob->soft; Mesh *me= ob->data; MFace *mface; BodyFace *bodyface; int a; /* alloc and copy faces*/ bodyface = sb->scratch->bodyface = MEM_mallocN(sizeof(BodyFace)*me->totface, "SB_body_Faces"); //memcpy(sb->scratch->mface, me->mface, sizeof(MFace)*me->totface); mface = me->mface; for (a=0; atotface; a++, mface++, bodyface++) { bodyface->v1 = mface->v1; bodyface->v2 = mface->v2; bodyface->v3 = mface->v3; bodyface->v4 = mface->v4; bodyface->ext_force[0] = bodyface->ext_force[1] = bodyface->ext_force[2] = 0.0f; bodyface->flag =0; } sb->scratch->totface = me->totface; } static void reference_to_scratch(Object *ob) { SoftBody *sb= ob->soft; ReferenceVert *rp; BodyPoint *bp; float accu_pos[3] ={0.f, 0.f, 0.f}; float accu_mass = 0.f; int a; sb->scratch->Ref.ivert = MEM_mallocN(sizeof(ReferenceVert)*sb->totpoint, "SB_Reference"); bp= ob->soft->bpoint; rp= sb->scratch->Ref.ivert; for (a=0; atotpoint; a++, rp++, bp++) { copy_v3_v3(rp->pos, bp->pos); add_v3_v3(accu_pos, bp->pos); accu_mass += _final_mass(ob, bp); } mul_v3_fl(accu_pos, 1.0f/accu_mass); copy_v3_v3(sb->scratch->Ref.com, accu_pos); /* printf("reference_to_scratch\n"); */ } /* * helper function to get proper spring length * when object is rescaled */ static float globallen(float *v1, float *v2, Object *ob) { float p1[3], p2[3]; copy_v3_v3(p1, v1); mul_m4_v3(ob->obmat, p1); copy_v3_v3(p2, v2); mul_m4_v3(ob->obmat, p2); return len_v3v3(p1, p2); } static void makelatticesprings(Lattice *lt, BodySpring *bs, int dostiff, Object *ob) { BPoint *bp=lt->def, *bpu; int u, v, w, dv, dw, bpc=0, bpuc; dv= lt->pntsu; dw= dv*lt->pntsv; for (w=0; wpntsw; w++) { for (v=0; vpntsv; v++) { for (u=0, bpuc=0, bpu=NULL; upntsu; u++, bp++, bpc++) { if (w) { bs->v1 = bpc; bs->v2 = bpc-dw; bs->springtype=SB_EDGE; bs->len= globallen((bp-dw)->vec, bp->vec, ob); bs++; } if (v) { bs->v1 = bpc; bs->v2 = bpc-dv; bs->springtype=SB_EDGE; bs->len= globallen((bp-dv)->vec, bp->vec, ob); bs++; } if (u) { bs->v1 = bpuc; bs->v2 = bpc; bs->springtype=SB_EDGE; bs->len= globallen((bpu)->vec, bp->vec, ob); bs++; } if (dostiff) { if (w) { if ( v && u ) { bs->v1 = bpc; bs->v2 = bpc-dw-dv-1; bs->springtype=SB_BEND; bs->len= globallen((bp-dw-dv-1)->vec, bp->vec, ob); bs++; } if ( (v < lt->pntsv-1) && (u != 0) ) { bs->v1 = bpc; bs->v2 = bpc-dw+dv-1; bs->springtype=SB_BEND; bs->len= globallen((bp-dw+dv-1)->vec, bp->vec, ob); bs++; } } if (w < lt->pntsw -1) { if ( v && u ) { bs->v1 = bpc; bs->v2 = bpc+dw-dv-1; bs->springtype=SB_BEND; bs->len= globallen((bp+dw-dv-1)->vec, bp->vec, ob); bs++; } if ( (v < lt->pntsv-1) && (u != 0) ) { bs->v1 = bpc; bs->v2 = bpc+dw+dv-1; bs->springtype=SB_BEND; bs->len= globallen((bp+dw+dv-1)->vec, bp->vec, ob); bs++; } } } bpu = bp; bpuc = bpc; } } } } /* makes totally fresh start situation */ static void lattice_to_softbody(Scene *scene, Object *ob) { Lattice *lt= ob->data; SoftBody *sb; int totvert, totspring = 0, a; BodyPoint *bp; BPoint *bpnt = lt->def; totvert= lt->pntsu*lt->pntsv*lt->pntsw; if (ob->softflag & OB_SB_EDGES) { totspring = ((lt->pntsu - 1) * lt->pntsv + (lt->pntsv - 1) * lt->pntsu) * lt->pntsw + lt->pntsu*lt->pntsv * (lt->pntsw - 1); if (ob->softflag & OB_SB_QUADS) { totspring += 4 * (lt->pntsu - 1) * (lt->pntsv -1) * (lt->pntsw - 1); } } /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ renew_softbody(scene, ob, totvert, totspring); sb= ob->soft; /* can be created in renew_softbody() */ bp = sb->bpoint; /* same code used as for mesh vertices */ for (a = 0; a < totvert; a++, bp++, bpnt++) { if ((ob->softflag & OB_SB_GOAL) && sb->vertgroup) { get_scalar_from_vertexgroup(ob, a, (short) (sb->vertgroup - 1), &bp->goal); } else { if (ob->softflag & OB_SB_GOAL) { bp->goal = sb->defgoal; } } } /* create some helper edges to enable SB lattice to be useful at all */ if (ob->softflag & OB_SB_EDGES) { makelatticesprings(lt, ob->soft->bspring, ob->softflag & OB_SB_QUADS, ob); build_bps_springlist(ob); /* link bps to springs */ } } /* makes totally fresh start situation */ static void curve_surf_to_softbody(Scene *scene, Object *ob) { Curve *cu= ob->data; SoftBody *sb; BodyPoint *bp; BodySpring *bs; Nurb *nu; BezTriple *bezt; BPoint *bpnt; int a, curindex=0; int totvert, totspring = 0, setgoal=0; totvert= BKE_nurbList_verts_count(&cu->nurb); if (ob->softflag & OB_SB_EDGES) { if (ob->type==OB_CURVE) { totspring = totvert - BLI_listbase_count(&cu->nurb); } } /* renew ends with ob->soft with points and edges, also checks & makes ob->soft */ renew_softbody(scene, ob, totvert, totspring); sb= ob->soft; /* can be created in renew_softbody() */ /* set vars now */ bp= sb->bpoint; bs= sb->bspring; /* weights from bpoints, same code used as for mesh vertices */ /* if ((ob->softflag & OB_SB_GOAL) && sb->vertgroup) 2.4x hack*/ /* new! take the weights from curve vertex anyhow */ if (ob->softflag & OB_SB_GOAL) setgoal= 1; for (nu= cu->nurb.first; nu; nu= nu->next) { if (nu->bezt) { /* bezier case ; this is nicly said naive; who ever wrote this part, it was not me (JOW) :) */ /* a: never ever make tangent handles (sub) and or (ob)ject to collision */ /* b: rather calculate them using some C2 (C2= continuous in second derivate -> no jump in bending ) condition */ /* not too hard to do, but needs some more code to care for; some one may want look at it JOW 2010/06/12*/ for (bezt=nu->bezt, a=0; apntsu; a++, bezt++, bp+=3, curindex+=3) { if (setgoal) { bp->goal= bezt->weight; /* all three triples */ (bp+1)->goal= bp->goal; (bp+2)->goal= bp->goal; /*do not collide handles */ (bp+1)->loc_flag |= SBF_OUTOFCOLLISION; (bp+2)->loc_flag |= SBF_OUTOFCOLLISION; } if (totspring) { if (a>0) { bs->v1= curindex-3; bs->v2= curindex; bs->springtype=SB_HANDLE; bs->len = globallen((bezt-1)->vec[0], bezt->vec[0], ob); bs++; } bs->v1= curindex; bs->v2= curindex+1; bs->springtype=SB_HANDLE; bs->len = globallen(bezt->vec[0], bezt->vec[1], ob); bs++; bs->v1= curindex+1; bs->v2= curindex+2; bs->springtype=SB_HANDLE; bs->len = globallen(bezt->vec[1], bezt->vec[2], ob); bs++; } } } else { for (bpnt=nu->bp, a=0; apntsu*nu->pntsv; a++, bpnt++, bp++, curindex++) { if (setgoal) { bp->goal= bpnt->weight; } if (totspring && a>0) { bs->v1= curindex-1; bs->v2= curindex; bs->springtype=SB_EDGE; bs->len= globallen( (bpnt-1)->vec, bpnt->vec, ob ); bs++; } } } } if (totspring) { build_bps_springlist(ob); /* link bps to springs */ if (ob->softflag & OB_SB_SELF) {calculate_collision_balls(ob);} } } /* copies softbody result back in object */ static void softbody_to_object(Object *ob, float (*vertexCos)[3], int numVerts, int local) { SoftBody *sb= ob->soft; if (sb) { BodyPoint *bp= sb->bpoint; int a; if (sb->solverflags & SBSO_ESTIMATEIPO) {SB_estimate_transform(ob, sb->lcom, sb->lrot, sb->lscale);} /* inverse matrix is not uptodate... */ invert_m4_m4(ob->imat, ob->obmat); for (a=0; apos); if (local==0) mul_m4_v3(ob->imat, vertexCos[a]); /* softbody is in global coords, baked optionally not */ } } } /* +++ ************ maintaining scratch *************** */ static void sb_new_scratch(SoftBody *sb) { if (!sb) return; sb->scratch = MEM_callocN(sizeof(SBScratch), "SBScratch"); sb->scratch->colliderhash = BLI_ghash_ptr_new("sb_new_scratch gh"); sb->scratch->bodyface = NULL; sb->scratch->totface = 0; sb->scratch->aabbmax[0]=sb->scratch->aabbmax[1]=sb->scratch->aabbmax[2] = 1.0e30f; sb->scratch->aabbmin[0]=sb->scratch->aabbmin[1]=sb->scratch->aabbmin[2] = -1.0e30f; sb->scratch->Ref.ivert = NULL; } /* --- ************ maintaining scratch *************** */ /* ************ Object level, exported functions *************** */ /* allocates and initializes general main data */ SoftBody *sbNew(Scene *scene) { SoftBody *sb; sb= MEM_callocN(sizeof(SoftBody), "softbody"); sb->mediafrict= 0.5f; sb->nodemass= 1.0f; sb->grav= 9.8f; sb->physics_speed= 1.0f; sb->rklimit= 0.1f; sb->goalspring= 0.5f; sb->goalfrict= 0.0f; sb->mingoal= 0.0f; sb->maxgoal= 1.0f; sb->defgoal= 0.7f; sb->inspring= 0.5f; sb->infrict= 0.5f; /*todo backward file compat should copy inspring to inpush while reading old files*/ sb->inpush = 0.5f; sb->interval= 10; sb->sfra= scene->r.sfra; sb->efra= scene->r.efra; sb->colball = 0.49f; sb->balldamp = 0.50f; sb->ballstiff= 1.0f; sb->sbc_mode = 1; sb->minloops = 10; sb->maxloops = 300; sb->choke = 3; sb_new_scratch(sb); /*todo backward file compat should set sb->shearstiff = 1.0f while reading old files*/ sb->shearstiff = 1.0f; sb->solverflags |= SBSO_OLDERR; sb->pointcache = BKE_ptcache_add(&sb->ptcaches); if (!sb->effector_weights) sb->effector_weights = BKE_add_effector_weights(NULL); sb->last_frame= MINFRAME-1; return sb; } /* frees all */ void sbFree(SoftBody *sb) { free_softbody_intern(sb); BKE_ptcache_free_list(&sb->ptcaches); sb->pointcache = NULL; if (sb->effector_weights) MEM_freeN(sb->effector_weights); MEM_freeN(sb); } void sbFreeSimulation(SoftBody *sb) { free_softbody_intern(sb); } /* makes totally fresh start situation */ void sbObjectToSoftbody(Object *ob) { //ob->softflag |= OB_SB_REDO; free_softbody_intern(ob->soft); } static int object_has_edges(Object *ob) { if (ob->type==OB_MESH) { return ((Mesh*) ob->data)->totedge; } else if (ob->type==OB_LATTICE) { return 1; } else { return 0; } } /* SB global visible functions */ void sbSetInterruptCallBack(int (*f)(void)) { SB_localInterruptCallBack = f; } static void softbody_update_positions(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) { BodyPoint *bp; int a; if (!sb || !sb->bpoint) return; for (a=0, bp=sb->bpoint; aorigS, bp->origE); /* copy the position of the goals at desired end time */ copy_v3_v3(bp->origE, vertexCos[a]); /* vertexCos came from local world, go global */ mul_m4_v3(ob->obmat, bp->origE); /* just to be save give bp->origT a defined value * will be calculated in interpolate_exciter() */ copy_v3_v3(bp->origT, bp->origE); } } /* void SB_estimate_transform */ /* input Object *ob out (says any object that can do SB like mesh, lattice, curve ) * output float lloc[3], float lrot[3][3], float lscale[3][3] * that is: * a precise position vector denoting the motion of the center of mass * give a rotation/scale matrix using averaging method, that's why estimate and not calculate * see: this is kind of reverse engineering: having to states of a point cloud and recover what happend * our advantage here we know the identity of the vertex * there are others methods giving other results. * lloc, lrot, lscale are allowed to be NULL, just in case you don't need it. * should be pretty useful for pythoneers :) * not! velocity .. 2nd order stuff * vcloud_estimate_transform see */ void SB_estimate_transform(Object *ob, float lloc[3], float lrot[3][3], float lscale[3][3]) { BodyPoint *bp; ReferenceVert *rp; SoftBody *sb = NULL; float (*opos)[3]; float (*rpos)[3]; float com[3], rcom[3]; int a; if (!ob ||!ob->soft) return; /* why did we get here ? */ sb= ob->soft; if (!sb || !sb->bpoint) return; opos= MEM_callocN((sb->totpoint)*3*sizeof(float), "SB_OPOS"); rpos= MEM_callocN((sb->totpoint)*3*sizeof(float), "SB_RPOS"); /* might filter vertex selection with a vertex group */ for (a=0, bp=sb->bpoint, rp=sb->scratch->Ref.ivert; atotpoint; a++, bp++, rp++) { copy_v3_v3(rpos[a], rp->pos); copy_v3_v3(opos[a], bp->pos); } vcloud_estimate_transform(sb->totpoint, opos, NULL, rpos, NULL, com, rcom, lrot, lscale); //sub_v3_v3(com, rcom); if (lloc) copy_v3_v3(lloc, com); copy_v3_v3(sb->lcom, com); if (lscale) copy_m3_m3(sb->lscale, lscale); if (lrot) copy_m3_m3(sb->lrot, lrot); MEM_freeN(opos); MEM_freeN(rpos); } static void softbody_reset(Object *ob, SoftBody *sb, float (*vertexCos)[3], int numVerts) { BodyPoint *bp; int a; for (a=0, bp=sb->bpoint; apos, vertexCos[a]); mul_m4_v3(ob->obmat, bp->pos); /* yep, sofbody is global coords*/ copy_v3_v3(bp->origS, bp->pos); copy_v3_v3(bp->origE, bp->pos); copy_v3_v3(bp->origT, bp->pos); bp->vec[0] = bp->vec[1] = bp->vec[2] = 0.0f; /* the bp->prev*'s are for rolling back from a canceled try to propagate in time * adaptive step size algo in a nutshell: * 1. set scheduled time step to new dtime * 2. try to advance the scheduled time step, being optimistic execute it * 3. check for success * 3.a we 're fine continue, may be we can increase scheduled time again ?? if so, do so! * 3.b we did exceed error limit --> roll back, shorten the scheduled time and try again at 2. * 4. check if we did reach dtime * 4.a nope we need to do some more at 2. * 4.b yup we're done */ copy_v3_v3(bp->prevpos, bp->pos); copy_v3_v3(bp->prevvec, bp->vec); copy_v3_v3(bp->prevdx, bp->vec); copy_v3_v3(bp->prevdv, bp->vec); } /* make a nice clean scratch struc */ free_scratch(sb); /* clear if any */ sb_new_scratch(sb); /* make a new */ sb->scratch->needstobuildcollider=1; zero_v3(sb->lcom); unit_m3(sb->lrot); unit_m3(sb->lscale); /* copy some info to scratch */ /* we only need that if we want to reconstruct IPO */ if (1) { reference_to_scratch(ob); SB_estimate_transform(ob, NULL, NULL, NULL); SB_estimate_transform(ob, NULL, NULL, NULL); } switch (ob->type) { case OB_MESH: if (ob->softflag & OB_SB_FACECOLL) mesh_faces_to_scratch(ob); break; case OB_LATTICE: break; case OB_CURVE: case OB_SURF: break; default: break; } } static void softbody_step(Scene *scene, Object *ob, SoftBody *sb, float dtime) { /* the simulator */ float forcetime; double sct, sst; sst=PIL_check_seconds_timer(); /* Integration back in time is possible in theory, but pretty useless here. * So we refuse to do so. Since we do not know anything about 'outside' changes * especially colliders we refuse to go more than 10 frames. */ if (dtime < 0 || dtime > 10.5f) return; ccd_update_deflector_hash(scene, ob, sb->scratch->colliderhash); if (sb->scratch->needstobuildcollider) { if (query_external_colliders(scene, ob)) { ccd_build_deflector_hash(scene, ob, sb->scratch->colliderhash); } sb->scratch->needstobuildcollider=0; } if (sb->solver_ID < 2) { /* special case of 2nd order Runge-Kutta type AKA Heun */ int mid_flags=0; float err = 0; float forcetimemax = 1.0f; /* set defaults guess we shall do one frame */ float forcetimemin = 0.01f; /* set defaults guess 1/100 is tight enough */ float timedone =0.0; /* how far did we get without violating error condition */ /* loops = counter for emergency brake * we don't want to lock up the system if physics fail */ int loops = 0; SoftHeunTol = sb->rklimit; /* humm .. this should be calculated from sb parameters and sizes */ /* adjust loop limits */ if (sb->minloops > 0) forcetimemax = dtime / sb->minloops; if (sb->maxloops > 0) forcetimemin = dtime / sb->maxloops; if (sb->solver_ID>0) mid_flags |= MID_PRESERVE; forcetime = forcetimemax; /* hope for integrating in one step */ while ( (ABS(timedone) < ABS(dtime)) && (loops < 2000) ) { /* set goals in time */ interpolate_exciter(ob, 200, (int)(200.0f*(timedone/dtime))); sb->scratch->flag &= ~SBF_DOFUZZY; /* do predictive euler step */ softbody_calc_forces(scene, ob, forcetime, timedone/dtime, 0); softbody_apply_forces(ob, forcetime, 1, NULL, mid_flags); /* crop new slope values to do averaged slope step */ softbody_calc_forces(scene, ob, forcetime, timedone/dtime, 0); softbody_apply_forces(ob, forcetime, 2, &err, mid_flags); softbody_apply_goalsnap(ob); if (err > SoftHeunTol) { /* error needs to be scaled to some quantity */ if (forcetime > forcetimemin) { forcetime = max_ff(forcetime / 2.0f, forcetimemin); softbody_restore_prev_step(ob); //printf("down, "); } else { timedone += forcetime; } } else { float newtime = forcetime * 1.1f; /* hope for 1.1 times better conditions in next step */ if (sb->scratch->flag & SBF_DOFUZZY) { //if (err > SoftHeunTol/(2.0f*sb->fuzzyness)) { /* stay with this stepsize unless err really small */ newtime = forcetime; //} } else { if (err > SoftHeunTol/2.0f) { /* stay with this stepsize unless err really small */ newtime = forcetime; } } timedone += forcetime; newtime = min_ff(forcetimemax, max_ff(newtime, forcetimemin)); //if (newtime > forcetime) printf("up, "); if (forcetime > 0.0f) forcetime = min_ff(dtime - timedone, newtime); else forcetime = max_ff(dtime - timedone, newtime); } loops++; if (sb->solverflags & SBSO_MONITOR ) { sct = PIL_check_seconds_timer(); if (sct - sst > 0.5) printf("%3.0f%% \r", 100.0f * timedone / dtime); } /* ask for user break */ if (SB_localInterruptCallBack && SB_localInterruptCallBack()) break; } /* move snapped to final position */ interpolate_exciter(ob, 2, 2); softbody_apply_goalsnap(ob); // if (G.debug & G_DEBUG) { if (sb->solverflags & SBSO_MONITOR ) { if (loops > HEUNWARNLIMIT) /* monitor high loop counts */ printf("\r needed %d steps/frame", loops); } } else if (sb->solver_ID == 2) { /* do semi "fake" implicit euler */ //removed }/*SOLVER SELECT*/ else if (sb->solver_ID == 4) { /* do semi "fake" implicit euler */ }/*SOLVER SELECT*/ else if (sb->solver_ID == 3) { /* do "stupid" semi "fake" implicit euler */ //removed }/*SOLVER SELECT*/ else { printf("softbody no valid solver ID!"); }/*SOLVER SELECT*/ if (sb->plastic) { apply_spring_memory(ob);} if (sb->solverflags & SBSO_MONITOR ) { sct=PIL_check_seconds_timer(); if ((sct - sst > 0.5) || (G.debug & G_DEBUG)) printf(" solver time %f sec %s\n", sct-sst, ob->id.name); } } /* simulates one step. framenr is in frames */ void sbObjectStep(Scene *scene, Object *ob, float cfra, float (*vertexCos)[3], int numVerts) { SoftBody *sb= ob->soft; PointCache *cache; PTCacheID pid; float dtime, timescale; int framedelta, framenr, startframe, endframe; int cache_result; cache= sb->pointcache; framenr= (int)cfra; framedelta= framenr - cache->simframe; BKE_ptcache_id_from_softbody(&pid, ob, sb); BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, ×cale); /* check for changes in mesh, should only happen in case the mesh * structure changes during an animation */ if (sb->bpoint && numVerts != sb->totpoint) { BKE_ptcache_invalidate(cache); return; } /* clamp frame ranges */ if (framenr < startframe) { BKE_ptcache_invalidate(cache); return; } else if (framenr > endframe) { framenr = endframe; } /* verify if we need to create the softbody data */ if (sb->bpoint == NULL || ((ob->softflag & OB_SB_EDGES) && !ob->soft->bspring && object_has_edges(ob))) { switch (ob->type) { case OB_MESH: mesh_to_softbody(scene, ob); break; case OB_LATTICE: lattice_to_softbody(scene, ob); break; case OB_CURVE: case OB_SURF: curve_surf_to_softbody(scene, ob); break; default: renew_softbody(scene, ob, numVerts, 0); break; } softbody_update_positions(ob, sb, vertexCos, numVerts); softbody_reset(ob, sb, vertexCos, numVerts); } /* still no points? go away */ if (sb->totpoint==0) { return; } if (framenr == startframe) { BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED); /* first frame, no simulation to do, just set the positions */ softbody_update_positions(ob, sb, vertexCos, numVerts); BKE_ptcache_validate(cache, framenr); cache->flag &= ~PTCACHE_REDO_NEEDED; sb->last_frame = framenr; return; } /* try to read from cache */ cache_result = BKE_ptcache_read(&pid, (float)framenr+scene->r.subframe); if (cache_result == PTCACHE_READ_EXACT || cache_result == PTCACHE_READ_INTERPOLATED) { softbody_to_object(ob, vertexCos, numVerts, sb->local); BKE_ptcache_validate(cache, framenr); if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED) BKE_ptcache_write(&pid, framenr); sb->last_frame = framenr; return; } else if (cache_result==PTCACHE_READ_OLD) { ; /* do nothing */ } else if (/*ob->id.lib || */(cache->flag & PTCACHE_BAKED)) { /* "library linking & pointcaches" has to be solved properly at some point */ /* if baked and nothing in cache, do nothing */ BKE_ptcache_invalidate(cache); return; } if (framenr!=sb->last_frame+1) return; /* if on second frame, write cache for first frame */ if (cache->simframe == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact==0)) BKE_ptcache_write(&pid, startframe); softbody_update_positions(ob, sb, vertexCos, numVerts); /* checking time: */ dtime = framedelta*timescale; /* do simulation */ softbody_step(scene, ob, sb, dtime); softbody_to_object(ob, vertexCos, numVerts, 0); BKE_ptcache_validate(cache, framenr); BKE_ptcache_write(&pid, framenr); sb->last_frame = framenr; }