archive/blender-archive
Archived
1
1
Fork 0
Code Pull Requests Packages Activity
This repository has been archived on 2023-10-09. You can view files and clone it, but cannot push or open issues or pull requests.
Files
6bf9c8ac8918a6656c20cca2c08f20b38e90980c
blender-archive/source/blender/blenkernel/intern/softbody.c
Dalai Felinto aeb8e81f27 Render Layers and Collections (merge from render-layers) 
Design Documents
----------------

* https://wiki.blender.org/index.php/Dev:2.8/Source/Layers

* https://wiki.blender.org/index.php/Dev:2.8/Source/DataDesignRevised

User Commit Log
---------------

* New Layer and Collection system to replace render layers and viewport layers.

* A layer is a set of collections of objects (and their drawing options) required for specific tasks.

* A collection is a set of objects, equivalent of the old layers in Blender. A collection can be shared across multiple layers.

* All Scenes have a master collection that all other collections are children of.

* New collection "context" tab (in Properties Editor)

* New temporary viewport "collections" panel to control per-collection
visibility

Missing User Features
---------------------

* Collection "Filter"
  Option to add objects based on their names

* Collection Manager operators
  The existing buttons  are placeholders

* Collection Manager drawing
  The editor main region is empty

* Collection Override

* Per-Collection engine settings
  This will come as a separate commit, as part of the clay-engine branch

Dev Commit Log
--------------

* New DNA file (DNA_layer_types.h) with the new structs
  We are replacing Base by a new extended Base while keeping it backward
  compatible with some legacy settings (i.e., lay, flag_legacy).

  Renamed all Base to BaseLegacy to make it clear the areas of code that
  still need to be converted

  Note: manual changes were required on - deg_builder_nodes.h, rna_object.c, KX_Light.cpp

* Unittesting for main syncronization requirements
  - read, write, add/copy/remove objects, copy scene, collection
  link/unlinking, context)

* New Editor: Collection Manager
  Based on patch by Julian Eisel
  This is extracted from the layer-manager branch. With the following changes:

    - Renamed references of layer manager to collections manager

    - I doesn't include the editors/space_collections/ draw and util files

    - The drawing code itself will be implemented separately by Julian

* Base / Object:
  A little note about them. Original Blender code would try to keep them
  in sync through the code, juggling flags back and forth. This will now
  be handled by Depsgraph, keeping Object and Bases more separated
  throughout the non-rendering code.

  Scene.base is being cleared in doversion, and the old viewport drawing
  code was poorly converted to use the new bases while the new viewport
  code doesn't get merged and replace the old one.

Python API Changes
------------------

```
- scene.layers
+ # no longer exists

- scene.objects
+ scene.scene_layers.active.objects

- scene.objects.active
+ scene.render_layers.active.objects.active

- bpy.context.scene.objects.link()
+ bpy.context.scene_collection.objects.link()

- bpy_extras.object_utils.object_data_add(context, obdata, operator=None, use_active_layer=True, name=None)
+ bpy_extras.object_utils.object_data_add(context, obdata, operator=None, name=None)

- bpy.context.object.select
+ bpy.context.object.select = True
+ bpy.context.object.select = False
+ bpy.context.object.select_get()
+ bpy.context.object.select_set(action='SELECT')
+ bpy.context.object.select_set(action='DESELECT')

-AddObjectHelper.layers
+ # no longer exists
```
2017-02-07 11:11:00 +01:00

3770 lines
109 KiB
C
Raw Blame History

/*
* ***** 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 <math.h>
#include <stdlib.h>
#include <string.h>
#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 "DNA_group_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"
/* 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 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 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 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 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 group overrides scene when not NULL.
*/
static void ccd_build_deflector_hash(Scene *scene, Group *group, Object *vertexowner, GHash *hash)
{
Object *ob;
if (!hash) return;
if (group) {
/* Explicit collision group */
for (GroupObject *go = group->gobject.first; go; go = go->next) {
ob = go->ob;
if (ob == vertexowner || ob->type != OB_MESH)
continue;
ccd_build_deflector_hash_single(hash, ob);
}
}
else {
for (BaseLegacy *base = scene->base.first; base; base = base->next) {
/*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 */
continue;
}
ccd_build_deflector_hash_single(hash, ob);
}
}
}
}
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 group overrides scene when not NULL.
*/
static void ccd_update_deflector_hash(Scene *scene, Group *group, Object *vertexowner, GHash *hash)
{
Object *ob;
if ((!hash) || (!vertexowner)) return;
if (group) {
/* Explicit collision group */
for (GroupObject *go = group->gobject.first; go; go = go->next) {
ob = go->ob;
if (ob == vertexowner || ob->type != OB_MESH)
continue;
ccd_update_deflector_hash_single(hash, ob);
}
}
else {
for (BaseLegacy *base = scene->base.first; base; base = base->next) {
/*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 */
continue;
}
ccd_update_deflector_hash_single(hash, ob);
}
}
}
}
/*--- 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 {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; i<totpoint; i++) {
BodyPoint *bp = &sb->bpoint[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; k<sb->totkey; 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
*/
/* +++ dependency information functions*/
/**
* \note group overrides scene when not NULL.
*/
static bool are_there_deflectors(Scene *scene, Group *group, unsigned int layer)
{
if (group) {
for (GroupObject *go = group->gobject.first; go; go = go->next) {
if (go->ob->pd && go->ob->pd->deflect)
return 1;
}
}
else {
for (BaseLegacy *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, Group *group, Object *me)
{
return(are_there_deflectors(scene, group, 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 (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*/
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) {
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*/
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);
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)<outerfacethickness) {
if (isect_point_tri_prism_v3(nv1, face_v1, face_v2, face_v3) ) {
float df;
if (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], unsigned int UNUSED(par_layer), 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*/
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))
{
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; a<sb->scratch->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->lay, 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->lay, 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; a<sb->scratch->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], unsigned int UNUSED(par_layer), 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*/
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))
{
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 *do_effector)
{
SoftBody *sb = ob->soft;
int a;
float damp;
float feedback[3];
if (sb && sb->totspring) {
for (a=ifirst; a<ilast; a++) {
BodySpring *bs = &sb->bspring[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);
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 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; i<totthread; i++) {
sb_threads[i].scene = scene;
sb_threads[i].ob = ob;
sb_threads[i].forcetime = 0.0; // not used here
sb_threads[i].timenow = timenow;
sb_threads[i].ilast = left;
left = left - dec;
if (left >0) {
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; i<totthread; i++)
BLI_insert_thread(&threads, &sb_threads[i]);
BLI_end_threads(&threads);
}
else
exec_scan_for_ext_spring_forces(&sb_threads[0]);
/* clean up */
MEM_freeN(sb_threads);
pdEndEffectors(&do_effector);
}
/* --- the spring external section*/
static int choose_winner(float*w, float* pos, float*a, float*b, float*c, float*ca, float*cb, float*cc)
{
float mindist, cp;
int winner =1;
mindist = fabsf(dot_v3v3(pos, a));
cp = fabsf(dot_v3v3(pos, b));
if ( mindist < cp ) {
mindist = cp;
winner =2;
}
cp = fabsf(dot_v3v3(pos, c));
if (mindist < cp ) {
mindist = cp;
winner =3;
}
switch (winner) {
case 1: copy_v3_v3(w, ca); break;
case 2: copy_v3_v3(w, cb); break;
case 3: copy_v3_v3(w, cc);
}
return(winner);
}
static int sb_detect_vertex_collisionCached(
float opco[3], float facenormal[3], float *damp,
float force[3], unsigned int UNUSED(par_layer), struct Object *vertexowner,
float time, float vel[3], float *intrusion)
{
Object *ob= NULL;
GHash *hash;
GHashIterator *ihash;
float nv1[3], nv2[3], nv3[3], edge1[3], edge2[3], d_nvect[3], dv1[3], ve[3], avel[3] = {0.0, 0.0, 0.0},
vv1[3], vv2[3], vv3[3], coledge[3] = {0.0f, 0.0f, 0.0f}, mindistedge = 1000.0f,
outerforceaccu[3], innerforceaccu[3],
facedist, /* n_mag, */ /* UNUSED */ force_mag_norm, minx, miny, minz, maxx, maxy, maxz,
innerfacethickness = -0.5f, outerfacethickness = 0.2f,
ee = 5.0f, ff = 0.1f, fa=1;
int a, deflected=0, cavel=0, ci=0;
/* init */
*intrusion = 0.0f;
hash = vertexowner->soft->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*/
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))
{
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->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]));
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);
}
*/
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 */
/**/
printf("bodypoint <bpi> 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;
}
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 *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 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);
/* exploit force(a, b) == -force(b, a) part2/2 */
sub_v3_v3v3(dvel, velcenter, obp->vec);
mul_v3_fl(dvel, _final_mass(ob, bp));
madd_v3_v3fl(obp->force, dvel, sb->balldamp);
madd_v3_v3fl(obp->force, def, -f * (1.0f - 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]);
/* 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);
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->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; i<totthread; i++) {
sb_threads[i].scene = scene;
sb_threads[i].ob = ob;
sb_threads[i].forcetime = forcetime;
sb_threads[i].timenow = timenow;
sb_threads[i].ilast = left;
left = left - dec;
if (left >0) {
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; i<totthread; i++)
BLI_insert_thread(&threads, &sb_threads[i]);
BLI_end_threads(&threads);
}
else
exec_softbody_calc_forces(&sb_threads[0]);
/* clean up */
MEM_freeN(sb_threads);
}
static void softbody_calc_forcesEx(Scene *scene, Object *ob, float forcetime, float timenow)
{
/* 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 *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, sb->collision_group, 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)
{
/* redirection to the new threaded Version */
if (!(G.debug_value & 0x10)) { // 16
softbody_calc_forcesEx(scene, ob, forcetime, timenow);
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 do_deflector, do_selfcollision, do_springcollision, do_aero;
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, sb->collision_group, 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);
}
bp = sb->bpoint;
for (int a = sb->totpoint; a > 0; a--, 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;
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 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 (( 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));
madd_v3_v3fl(bp->force, def, f * (1.0f - sb->balldamp));
madd_v3_v3fl(bp->force, dvel, sb->balldamp);
/* 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);
madd_v3_v3fl(obp->force, dvel, sb->balldamp);
madd_v3_v3fl(obp->force, def, -f * (1.0f - 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 */
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 */
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 */
}
/* +++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 (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
*/
madd_v3_v3fl(bp->pos, facenormal, -intrusion);
}
else {
sub_v3_v3v3(cfforce, bp->vec, vel);
madd_v3_v3fl(bp->force, cfforce, -cf * 50.0f);
}
sb->scratch->flag |= SBF_DOFUZZY;
bp->loc_flag |= SBF_DOFUZZY;
bp->choke = sb->choke*0.01f;
}
else {
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 */
if (ob->softflag & OB_SB_EDGES) {
if (sb->bspring) { /* spring list exists at all ? */
for (int 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)
// rather remove nl_falgs from code .. will make things a lot cleaner
sb_spring_force(ob, sb->totpoint-a, bs, iks, forcetime);
}/* 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);
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)<freezeloc )&&(dot_v3v3(bp->force, bp->force)<freezeforce )) {
bp->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 <viscosity>
* 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; a<sb->totspring; 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 calcualtions here
* will be overwritten sbObjectStep() to receive
* actual modifier stack positions
*/
if (me->totvert) {
bp= ob->soft->bpoint;
for (a=0; a<me->totvert; 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; a<sb->totspring; 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; a<me->totvert; 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; a<sb->totpoint; 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; w<lt->pntsw; w++) {
for (v=0; v<lt->pntsv; v++) {
for (u=0, bpuc=0, bpu=NULL; u<lt->pntsu; 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(&lt->dvert[a], defgroup_index);
}
else {
bp->goal *= bpnt->weight;
}
if (defgroup_index_mass != -1) {
bp->mass *= defvert_find_weight(&lt->dvert[a], defgroup_index_mass);
}
if (defgroup_index_spring != -1) {
bp->springweight *= defvert_find_weight(&lt->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; a<nu->pntsu; 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; a<nu->pntsu*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; a<numVerts; a++, bp++) {
copy_v3_v3(vertexCos[a], bp->pos);
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; a<numVerts; a++, bp++) {
/* store where goals are now */
copy_v3_v3(bp->origS, 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; a<sb->totpoint; 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; a<numVerts; a++, bp++) {
copy_v3_v3(bp->pos, 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, sb->collision_group, ob, sb->scratch->colliderhash);
if (sb->scratch->needstobuildcollider) {
if (query_external_colliders(scene, sb->collision_group, ob)) {
ccd_build_deflector_hash(scene, 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(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(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 {
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, &timescale);
/* 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 */
bool can_simulate = (framenr == sb->last_frame + 1) && !(cache->flag & PTCACHE_BAKED);
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)
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 (!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(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;
}
View Git Blame Copy Permalink