/*
* 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.
*/
/** \file
* \ingroup bke
*/
/**
* variables on the UI for now
*
* float mediafrict; friction to env
* float nodemass; softbody mass of *vertex*
* float grav; softbody amount of gravitation 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 "CLG_log.h"
#include "MEM_guardedalloc.h"
/* types */
#include "DNA_collection_types.h"
#include "DNA_curve_types.h"
#include "DNA_lattice_types.h"
#include "DNA_meshdata_types.h"
#include "DNA_mesh_types.h"
#include "DNA_object_types.h"
#include "DNA_scene_types.h"
#include "BLI_math.h"
#include "BLI_utildefines.h"
#include "BLI_listbase.h"
#include "BLI_ghash.h"
#include "BLI_threads.h"
#include "BKE_collection.h"
#include "BKE_collision.h"
#include "BKE_curve.h"
#include "BKE_effect.h"
#include "BKE_global.h"
#include "BKE_layer.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 "DEG_depsgraph.h"
#include "DEG_depsgraph_query.h"
#include "PIL_time.h"
static CLG_LogRef LOG = {"bke.softbody"};
/* 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;
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 initial 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 *effectors;
int do_deflector;
float fieldfactor;
float windfactor;
int nr;
int tot;
} SB_thread_context;
#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);
/*+++ 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 happened 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);
}
}
CLOG_ERROR(&LOG, "sb or bp == NULL");
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);
}
}
CLOG_ERROR(&LOG, "sb or bp == NULL");
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 mvert_num, tri_num;
const MVert *mvert;
const MVert *mprevvert;
const MVertTri *tri;
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;
const MVertTri *vt;
float hull;
int i;
cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision);
/* first some paranoia checks */
if (!cmd) return NULL;
if (!cmd->mvert_num || !cmd->tri_num) return NULL;
pccd_M = MEM_mallocN(sizeof(ccd_Mesh), "ccd_Mesh");
pccd_M->mvert_num = cmd->mvert_num;
pccd_M->tri_num = cmd->tri_num;
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->mvert_num; i++) {
const float *v;
/* evaluate limits */
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->tri = MEM_dupallocN(cmd->tri);
/* OBBs for idea1 */
pccd_M->mima = MEM_mallocN(sizeof(ccdf_minmax) * pccd_M->tri_num, "ccd_Mesh_Faces_mima");
/* anyhoo we need to walk the list of faces and find OBB they live in */
for (i = 0, mima = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
const float *v;
mima->minx=mima->miny=mima->minz=1e30f;
mima->maxx=mima->maxy=mima->maxz=-1e30f;
v = pccd_M->mvert[vt->tri[0]].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);
v = pccd_M->mvert[vt->tri[1]].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);
v = pccd_M->mvert[vt->tri[2]].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);
}
return pccd_M;
}
static void ccd_mesh_update(Object *ob, ccd_Mesh *pccd_M)
{
CollisionModifierData *cmd;
ccdf_minmax *mima;
const MVertTri *vt;
float hull;
int i;
cmd =(CollisionModifierData *)modifiers_findByType(ob, eModifierType_Collision);
/* first some paranoia checks */
if (!cmd) return;
if (!cmd->mvert_num || !cmd->tri_num) return;
if ((pccd_M->mvert_num != cmd->mvert_num) ||
(pccd_M->tri_num != cmd->tri_num))
{
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((void *)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->mvert_num; i++) {
const float *v;
/* evaluate limits */
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 */
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);
}
/* anyhoo we need to walk the list of faces and find OBB they live in */
for (i = 0, mima = pccd_M->mima, vt = pccd_M->tri; i < pccd_M->tri_num; i++, mima++, vt++) {
const float *v;
mima->minx=mima->miny=mima->minz=1e30f;
mima->maxx=mima->maxy=mima->maxz=-1e30f;
/* mvert */
v = pccd_M->mvert[vt->tri[0]].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);
v = pccd_M->mvert[vt->tri[1]].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);
v = pccd_M->mvert[vt->tri[2]].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);
/* mprevvert */
v = pccd_M->mprevvert[vt->tri[0]].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);
v = pccd_M->mprevvert[vt->tri[1]].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);
v = pccd_M->mprevvert[vt->tri[2]].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);
}
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((void *)ccdm->mvert);
MEM_freeN((void *)ccdm->tri);
if (ccdm->mprevvert) MEM_freeN((void *)ccdm->mprevvert);
MEM_freeN(ccdm->mima);
MEM_freeN(ccdm);
ccdm = NULL;
}
}
static void ccd_build_deflector_hash_single(GHash *hash, Object *ob)
{
/* only with deflecting set */
if (ob->pd && ob->pd->deflect) {
void **val_p;
if (!BLI_ghash_ensure_p(hash, ob, &val_p)) {
ccd_Mesh *ccdmesh = ccd_mesh_make(ob);
*val_p = ccdmesh;
}
}
}
/**
* \note collection overrides scene when not NULL.
*/
static void ccd_build_deflector_hash(Depsgraph *depsgraph, Collection *collection, Object *vertexowner, GHash *hash)
{
if (!hash) return;
unsigned int numobjects;
Object **objects = BKE_collision_objects_create(depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
for (int i = 0; i < numobjects; i++) {
Object *ob = objects[i];
if (ob->type == OB_MESH) {
ccd_build_deflector_hash_single(hash, ob);
}
}
BKE_collision_objects_free(objects);
}
static void ccd_update_deflector_hash_single(GHash *hash, Object *ob)
{
if (ob->pd && ob->pd->deflect) {
ccd_Mesh *ccdmesh = BLI_ghash_lookup(hash, ob);
if (ccdmesh) {
ccd_mesh_update(ob, ccdmesh);
}
}
}
/**
* \note collection overrides scene when not NULL.
*/
static void ccd_update_deflector_hash(Depsgraph *depsgraph, Collection *collection, Object *vertexowner, GHash *hash)
{
if ((!hash) || (!vertexowner)) return;
unsigned int numobjects;
Object **objects = BKE_collision_objects_create(depsgraph, vertexowner, collection, &numobjects, eModifierType_Collision);
for (int i = 0; i < numobjects; i++) {
Object *ob = objects[i];
if (ob->type == OB_MESH) {
ccd_update_deflector_hash_single(hash, ob);
}
}
BKE_collision_objects_free(objects);
}
/*--- collider caching and dicing ---*/
static int count_mesh_quads(Mesh *me)
{
int a, result = 0;
const MPoly *mp = me->mpoly;
if (mp) {
for (a = me->totpoly; a > 0; a--, mp++) {
if (mp->totloop == 4) {
result++;
}
}
}
return result;
}
static void add_mesh_quad_diag_springs(Object *ob)
{
Mesh *me= ob->data;
/*BodyPoint *bp;*/ /*UNUSED*/
int a;
if (ob->soft) {
int nofquads;
//float s_shear = ob->soft->shearstiff*ob->soft->shearstiff;
nofquads = count_mesh_quads(me);
if (nofquads) {
const MLoop *mloop = me->mloop;
const MPoly *mp = me->mpoly;
BodySpring *bs;
/* resize spring-array to hold additional quad springs */
ob->soft->bspring = MEM_recallocN(ob->soft->bspring, sizeof(BodySpring) * (ob->soft->totspring + nofquads * 2));
/* fill the tail */
a = 0;
bs = &ob->soft->bspring[ob->soft->totspring];
/*bp= ob->soft->bpoint; */ /*UNUSED*/
for (a = me->totpoly; a > 0; a--, mp++) {
if (mp->totloop == 4) {
bs->v1 = mloop[mp->loopstart + 0].v;
bs->v2 = mloop[mp->loopstart + 2].v;
bs->springtype = SB_STIFFQUAD;
bs++;
bs->v1 = mloop[mp->loopstart + 1].v;
bs->v2 = mloop[mp->loopstart + 3].v;
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 {
CLOG_ERROR(&LOG, "oops we should not get here");
}
}
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 strength 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
*/
/* +++ dependency information functions*/
/**
* \note collection overrides scene when not NULL.
*/
static int query_external_colliders(Depsgraph *depsgraph, Collection *collection)
{
unsigned int numobjects;
Object **objects = BKE_collision_objects_create(depsgraph, NULL, collection, &numobjects, eModifierType_Collision);
BKE_collision_objects_free(objects);
return (numobjects != 0);
}
/* --- dependency information functions*/
/* +++ the aabb "force" section*/
static int sb_detect_aabb_collisionCached(float UNUSED(force[3]), 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 (ccdm) {
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*/
CLOG_ERROR(&LOG, "missing cache error");
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], 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) {
const MVert *mvert= NULL;
const MVert *mprevvert= NULL;
if (ccdm) {
mvert = ccdm->mvert;
a = ccdm->mvert_num;
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*/
CLOG_ERROR(&LOG, "missing cache error");
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);
madd_v3_v3fl(nv1, mprevvert[a - 1].co, 1.0f - time);
}
/* 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);
madd_v3_v3fl(force, d_nvect, -df);
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], struct Object *vertexowner, float time)
{
Object *ob;
GHash *hash;
GHashIterator *ihash;
float nv1[3], nv2[3], nv3[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) {
const MVert *mvert = NULL;
const MVert *mprevvert = NULL;
const MVertTri *vt = NULL;
const ccdf_minmax *mima = NULL;
if (ccdm) {
mvert = ccdm->mvert;
vt = ccdm->tri;
mprevvert = ccdm->mprevvert;
mima = ccdm->mima;
a = ccdm->tri_num;
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*/
CLOG_ERROR(&LOG, "missing cache error");
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))
{
mima++;
vt++;
continue;
}
if (mvert) {
copy_v3_v3(nv1, mvert[vt->tri[0]].co);
copy_v3_v3(nv2, mvert[vt->tri[1]].co);
copy_v3_v3(nv3, mvert[vt->tri[2]].co);
if (mprevvert) {
mul_v3_fl(nv1, time);
madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
mul_v3_fl(nv2, time);
madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
mul_v3_fl(nv3, time);
madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
}
}
/* 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_segment_tri_v3(nv1, nv2, face_v1, face_v2, face_v3, &t, NULL) ||
isect_line_segment_tri_v3(nv2, nv3, face_v1, face_v2, face_v3, &t, NULL) ||
isect_line_segment_tri_v3(nv3, nv1, face_v1, face_v2, face_v3, &t, NULL) )
{
madd_v3_v3fl(force, d_nvect, -0.5f);
*damp=tune*ob->pd->pdef_sbdamp;
deflected = 2;
}
mima++;
vt++;
}/* 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;
zero_v3(feedback);
if (sb_detect_face_collisionCached(
sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos,
&damp, feedback, ob, timenow))
{
madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
// madd_v3_v3fl(bf->ext_force, feedback, tune);
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;
zero_v3(feedback);
if (sb_detect_face_pointCached(
sb->bpoint[bf->v1].pos, sb->bpoint[bf->v2].pos, sb->bpoint[bf->v3].pos,
&damp, feedback, ob, timenow))
{
madd_v3_v3fl(sb->bpoint[bf->v1].force, feedback, tune);
madd_v3_v3fl(sb->bpoint[bf->v2].force, feedback, tune);
madd_v3_v3fl(sb->bpoint[bf->v3].force, feedback, tune);
// madd_v3_v3fl(bf->ext_force, feedback, tune);
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);
}
}
}
}
/* --- 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], struct Object *vertexowner, float time)
{
Object *ob;
GHash *hash;
GHashIterator *ihash;
float nv1[3], nv2[3], nv3[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) {
const MVert *mvert = NULL;
const MVert *mprevvert = NULL;
const MVertTri *vt = NULL;
const ccdf_minmax *mima = NULL;
if (ccdm) {
mvert = ccdm->mvert;
mprevvert = ccdm->mprevvert;
vt = ccdm->tri;
mima = ccdm->mima;
a = ccdm->tri_num;
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*/
CLOG_ERROR(&LOG, "missing cache error");
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))
{
mima++;
vt++;
continue;
}
if (mvert) {
copy_v3_v3(nv1, mvert[vt->tri[0]].co);
copy_v3_v3(nv2, mvert[vt->tri[1]].co);
copy_v3_v3(nv3, mvert[vt->tri[2]].co);
if (mprevvert) {
mul_v3_fl(nv1, time);
madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
mul_v3_fl(nv2, time);
madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
mul_v3_fl(nv3, time);
madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
}
}
/* 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_segment_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;
madd_v3_v3fl(force, d_nvect, intrusiondepth);
*damp=ob->pd->pdef_sbdamp;
deflected = 2;
}
mima++;
vt++;
}/* 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 *effectors)
{
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, 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 (effectors) {
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);
BKE_effectors_apply(effectors, 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);
madd_v3_v3fl(bs->ext_force, vel, f * (1.0f - fabsf(dot_v3v3(vel, sp))));
}
else {
madd_v3_v3fl(bs->ext_force, vel, f); // to keep compatible with 2.45 release files
}
}
/* --- springs seeing wind */
}
}
}
}
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->effectors);
return NULL;
}
static void sb_sfesf_threads_run(struct Depsgraph *depsgraph, Scene *scene, struct Object *ob, float timenow, int totsprings, int *UNUSED(ptr_to_break_func(void)))
{
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 */
ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, ob->soft->effector_weights);
/* 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].effectors = effectors;
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_threadpool_init(&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) {
const MVert *mvert = NULL;
const MVert *mprevvert = NULL;
const MVertTri *vt = NULL;
const ccdf_minmax *mima = NULL;
if (ccdm) {
mvert = ccdm->mvert;
mprevvert = ccdm->mprevvert;
vt = ccdm->tri;
mima = ccdm->mima;
a = ccdm->tri_num;
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*/
CLOG_ERROR(&LOG, "missing cache error");
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))
{
mima++;
vt++;
continue;
}
if (mvert) {
copy_v3_v3(nv1, mvert[vt->tri[0]].co);
copy_v3_v3(nv2, mvert[vt->tri[1]].co);
copy_v3_v3(nv3, mvert[vt->tri[2]].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[vt->tri[0]].co);
sub_v3_v3v3(vv2, nv2, mprevvert[vt->tri[1]].co);
sub_v3_v3v3(vv3, nv3, mprevvert[vt->tri[2]].co);
mul_v3_fl(nv1, time);
madd_v3_v3fl(nv1, mprevvert[vt->tri[0]].co, 1.0f - time);
mul_v3_fl(nv2, time);
madd_v3_v3fl(nv2, mprevvert[vt->tri[1]].co, 1.0f - time);
mul_v3_fl(nv3, time);
madd_v3_v3fl(nv3, mprevvert[vt->tri[2]].co, 1.0f - time);
}
}
/* 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);
madd_v3_v3fl(outerforceaccu, d_nvect, force_mag_norm);
deflected = 3;
}
else {
madd_v3_v3fl(innerforceaccu, d_nvect, force_mag_norm);
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++;
}
}
mima++;
vt++;
}/* 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);
madd_v3_v3fl(force, coledge, force_mag_norm);
*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, time, vel, intrusion);
//deflected= sb_detect_vertex_collisionCachedEx(s_actpos, facenormal, cf, force, 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 */
#if 0
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]));
EIG_linear_solver_matrix_add(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];
EIG_linear_solver_matrix_add(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++) EIG_linear_solver_matrix_add(ia+i, op+ic+i, factor);
}
static void dfdv_goal(int ia, int ic, float factor)
{
int i;
for (i=0;i<3;i++) EIG_linear_solver_matrix_add(ia+i, ic+i, factor);
}
#endif /* if 0 */
static void sb_spring_force(Object *ob, int bpi, BodySpring *bs, float iks, float UNUSED(forcetime))
{
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 */
/**/
CLOG_WARN(&LOG, "bodypoint is not attached to spring <*bs>");
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;
}
madd_v3_v3fl(bp1->force, dir, (bs->len - distance) * forcefactor);
/* 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;
madd_v3_v3fl(bp1->force, dir, -kd);
}
/* 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 *effectors, 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;
/* initialize */
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 {
CLOG_ERROR(&LOG, "expected a SB here");
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;
/* running in a slice we must not assume anything done with obp neither alter the data of obp */
for (c=sb->totpoint, obp= sb->bpoint; c>0; c--, obp++) {
compare = (obp->colball + bp->colball);
sub_v3_v3v3(def, bp->pos, obp->pos);
/* rather check the AABBoxes before ever calculating the real distance */
/* mathematically it is completely 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));
madd_v3_v3fl(bp->force, def, f * (1.0f - sb->balldamp));
madd_v3_v3fl(bp->force, dvel, sb->balldamp);
}
}
}
}
/* 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]);
/* calculate 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 (effectors) {
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);
BKE_effectors_apply(effectors, 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);
madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
madd_v3_v3fl(bp->force, defforce, kd);
}
}
/* ---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)
sb_spring_force(ob, ilast-bb, bs, iks, forcetime);
}/* 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->effectors, 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 *effectors, 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].effectors = effectors;
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_threadpool_init(&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 */
/* 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(depsgraph, sb->collision_group);
/* 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(depsgraph, scene, ob, timenow, sb->totspring, NULL);
/* after spring scan because it uses Effoctors too */
ListBase *effectors = BKE_effectors_create(depsgraph, ob, NULL, sb->effector_weights);
if (do_deflector) {
float defforce[3];
do_deflector = sb_detect_aabb_collisionCached(defforce, ob, timenow);
}
sb_cf_threads_run(scene, ob, forcetime, timenow, sb->totpoint, NULL, effectors, 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 */
BKE_effectors_free(effectors);
}
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 0
if (sb->nodemass > 0.009999f) timeovermass = forcetime / sb->nodemass;
else timeovermass = forcetime / 0.009999f;
#endif
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
*/
/* 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 calculations 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;
int defgroup_index, defgroup_index_mass, defgroup_index_spring;
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;
defgroup_index = me->dvert ? (sb->vertgroup - 1) : -1;
defgroup_index_mass = me->dvert ? defgroup_name_index(ob, sb->namedVG_Mass) : -1;
defgroup_index_spring = me->dvert ? defgroup_name_index(ob, sb->namedVG_Spring_K) : -1;
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) {
BLI_assert(bp->goal == sb->defgoal);
}
if ((ob->softflag & OB_SB_GOAL) && (defgroup_index != -1)) {
bp->goal *= defvert_find_weight(&me->dvert[a], defgroup_index);
}
/* to proof the concept
* this enables per vertex *mass painting*
*/
if (defgroup_index_mass != -1) {
bp->mass *= defvert_find_weight(&me->dvert[a], defgroup_index_mass);
}
if (defgroup_index_spring != -1) {
bp->springweight *= defvert_find_weight(&me->dvert[a], defgroup_index_spring);
}
}
/* 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;
const Mesh *me = ob->data;
MLoopTri *looptri, *lt;
BodyFace *bodyface;
int a;
/* alloc and copy faces*/
sb->scratch->totface = poly_to_tri_count(me->totpoly, me->totloop);
looptri = lt = MEM_mallocN(sizeof(*looptri) * sb->scratch->totface, __func__);
BKE_mesh_recalc_looptri(
me->mloop, me->mpoly,
me->mvert,
me->totloop, me->totpoly,
looptri);
bodyface = sb->scratch->bodyface = MEM_mallocN(sizeof(BodyFace) * sb->scratch->totface, "SB_body_Faces");
for (a = 0; a < sb->scratch->totface; a++, lt++, bodyface++) {
bodyface->v1 = me->mloop[lt->tri[0]].v;
bodyface->v2 = me->mloop[lt->tri[1]].v;
bodyface->v3 = me->mloop[lt->tri[2]].v;
zero_v3(bodyface->ext_force);
bodyface->ext_force[0] = bodyface->ext_force[1] = bodyface->ext_force[2] = 0.0f;
bodyface->flag = 0;
}
MEM_freeN(looptri);
}
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;
int defgroup_index, defgroup_index_mass, defgroup_index_spring;
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;
defgroup_index = lt->dvert ? (sb->vertgroup - 1) : -1;
defgroup_index_mass = lt->dvert ? defgroup_name_index(ob, sb->namedVG_Mass) : -1;
defgroup_index_spring = lt->dvert ? defgroup_name_index(ob, sb->namedVG_Spring_K) : -1;
/* same code used as for mesh vertices */
for (a = 0; a < totvert; a++, bp++, bpnt++) {
if (ob->softflag & OB_SB_GOAL) {
BLI_assert(bp->goal == sb->defgoal);
}
if ((ob->softflag & OB_SB_GOAL) && (defgroup_index != -1)) {
bp->goal *= defvert_find_weight(<->dvert[a], defgroup_index);
}
else {
bp->goal *= bpnt->weight;
}
if (defgroup_index_mass != -1) {
bp->mass *= defvert_find_weight(<->dvert[a], defgroup_index_mass);
}
if (defgroup_index_spring != -1) {
bp->springweight *= defvert_find_weight(<->dvert[a], defgroup_index_spring);
}
}
/* 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->shared = MEM_callocN(sizeof(*sb->shared), "SoftBody_Shared");
sb->shared->pointcache = BKE_ptcache_add(&sb->shared->ptcaches);
if (!sb->effector_weights)
sb->effector_weights = BKE_effector_add_weights(NULL);
sb->last_frame= MINFRAME-1;
return sb;
}
/* frees all */
void sbFree(Object *ob)
{
SoftBody *sb = ob->soft;
if (sb == NULL) {
return;
}
free_softbody_intern(sb);
if ((ob->id.tag & LIB_TAG_NO_MAIN) == 0) {
/* Only free shared data on non-CoW copies */
BKE_ptcache_free_list(&sb->shared->ptcaches);
sb->shared->pointcache = NULL;
MEM_freeN(sb->shared);
}
if (sb->effector_weights)
MEM_freeN(sb->effector_weights);
MEM_freeN(sb);
ob->soft = NULL;
}
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 happened
* 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_v3 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_v3(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(struct Depsgraph *depsgraph, 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(depsgraph, sb->collision_group, ob, sb->scratch->colliderhash);
if (sb->scratch->needstobuildcollider) {
ccd_build_deflector_hash(depsgraph, sb->collision_group, 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(depsgraph, scene, ob, forcetime, timedone/dtime);
softbody_apply_forces(ob, forcetime, 1, NULL, mid_flags);
/* crop new slope values to do averaged slope step */
softbody_calc_forces(depsgraph, scene, ob, forcetime, timedone/dtime);
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 {
CLOG_ERROR(&LOG, "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);
}
}
static void sbStoreLastFrame(struct Depsgraph *depsgraph, Object *object, float framenr)
{
if (!DEG_is_active(depsgraph)) {
return;
}
Object *object_orig = DEG_get_original_object(object);
object->soft->last_frame = framenr;
object_orig->soft->last_frame = framenr;
}
/* simulates one step. framenr is in frames */
void sbObjectStep(struct Depsgraph *depsgraph, 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->shared->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;
sbStoreLastFrame(depsgraph, ob, framenr);
return;
}
/* try to read from cache */
bool can_write_cache = DEG_is_active(depsgraph);
bool can_simulate = (framenr == sb->last_frame + 1) && !(cache->flag & PTCACHE_BAKED) && can_write_cache;
cache_result = BKE_ptcache_read(&pid, (float)framenr+scene->r.subframe, can_simulate);
if (cache_result == PTCACHE_READ_EXACT || cache_result == PTCACHE_READ_INTERPOLATED ||
(!can_simulate && cache_result == PTCACHE_READ_OLD))
{
softbody_to_object(ob, vertexCos, numVerts, sb->local);
BKE_ptcache_validate(cache, framenr);
if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED && can_write_cache)
BKE_ptcache_write(&pid, framenr);
sbStoreLastFrame(depsgraph, ob, 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 */
if (can_write_cache) {
BKE_ptcache_invalidate(cache);
}
return;
}
if (!can_simulate)
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(depsgraph, scene, ob, sb, dtime);
softbody_to_object(ob, vertexCos, numVerts, 0);
BKE_ptcache_validate(cache, framenr);
BKE_ptcache_write(&pid, framenr);
sbStoreLastFrame(depsgraph, ob, framenr);
}