To have the 50% faster nearest_surface point.
Changed mesh_faces_nearest_point to return the face normal instead of collision normal
This commit is contained in:
Andre Susano Pinto
2008-08-13 19:22:35 +00:00
parent 6a8236a8da
commit 43bf03580f
8 changed files with 1046 additions and 457 deletions

View File

@@ -48,9 +48,6 @@
/* Math stuff for ray casting on mesh faces and for nearest surface */
static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest);
#define ISECT_EPSILON 1e-6
static float ray_tri_intersection(const BVHTreeRay *ray, const float m_dist, const float *v0, const float *v1, const float *v2)
{
float dist;
@@ -79,170 +76,324 @@ static float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, con
return FLT_MAX;
}
/*
* This calculates the distance from point to the plane
* Distance is negative if point is on the back side of plane
* Function adapted from David Eberly's distance tools (LGPL)
* http://www.geometrictools.com/LibFoundation/Distance/Distance.html
*/
static float point_plane_distance(const float *point, const float *plane_point, const float *plane_normal)
static float nearest_point_in_tri_surface(const float *v0,const float *v1,const float *v2,const float *p, int *v, int *e, float *nearest )
{
float pp[3];
VECSUB(pp, point, plane_point);
return INPR(pp, plane_normal);
}
static float choose_nearest(const float v0[2], const float v1[2], const float point[2], float closest[2])
{
float d[2][2], sdist[2];
VECSUB2D(d[0], v0, point);
VECSUB2D(d[1], v1, point);
float diff[3];
float e0[3];
float e1[3];
float A00;
float A01;
float A11;
float B0;
float B1;
float C;
float Det;
float S;
float T;
float sqrDist;
int lv = -1, le = -1;
VECSUB(diff, v0, p);
VECSUB(e0, v1, v0);
VECSUB(e1, v2, v0);
A00 = INPR ( e0, e0 );
A01 = INPR( e0, e1 );
A11 = INPR ( e1, e1 );
B0 = INPR( diff, e0 );
B1 = INPR( diff, e1 );
C = INPR( diff, diff );
Det = fabs( A00 * A11 - A01 * A01 );
S = A01 * B1 - A11 * B0;
T = A01 * B0 - A00 * B1;
sdist[0] = d[0][0]*d[0][0] + d[0][1]*d[0][1];
sdist[1] = d[1][0]*d[1][0] + d[1][1]*d[1][1];
if(sdist[0] < sdist[1])
if ( S + T <= Det )
{
if(closest)
VECCOPY2D(closest, v0);
return sdist[0];
if ( S < 0.0f )
{
if ( T < 0.0f ) // Region 4
{
if ( B0 < 0.0f )
{
T = 0.0f;
if ( -B0 >= A00 )
{
S = (float)1.0;
sqrDist = A00 + 2.0f * B0 + C;
lv = 1;
}
else
{
if(fabs(A00) > FLT_EPSILON)
S = -B0/A00;
else
S = 0.0f;
sqrDist = B0 * S + C;
le = 0;
}
}
else
{
S = 0.0f;
if ( B1 >= 0.0f )
{
T = 0.0f;
sqrDist = C;
lv = 0;
}
else if ( -B1 >= A11 )
{
T = 1.0f;
sqrDist = A11 + 2.0f * B1 + C;
lv = 2;
}
else
{
if(fabs(A11) > FLT_EPSILON)
T = -B1 / A11;
else
T = 0.0f;
sqrDist = B1 * T + C;
le = 1;
}
}
}
else // Region 3
{
S = 0.0f;
if ( B1 >= 0.0f )
{
T = 0.0f;
sqrDist = C;
lv = 0;
}
else if ( -B1 >= A11 )
{
T = 1.0f;
sqrDist = A11 + 2.0f * B1 + C;
lv = 2;
}
else
{
if(fabs(A11) > FLT_EPSILON)
T = -B1 / A11;
else
T = 0.0;
sqrDist = B1 * T + C;
le = 1;
}
}
}
else if ( T < 0.0f ) // Region 5
{
T = 0.0f;
if ( B0 >= 0.0f )
{
S = 0.0f;
sqrDist = C;
lv = 0;
}
else if ( -B0 >= A00 )
{
S = 1.0f;
sqrDist = A00 + 2.0f * B0 + C;
lv = 1;
}
else
{
if(fabs(A00) > FLT_EPSILON)
S = -B0 / A00;
else
S = 0.0f;
sqrDist = B0 * S + C;
le = 0;
}
}
else // Region 0
{
// Minimum at interior lv
float invDet;
if(fabs(Det) > FLT_EPSILON)
invDet = 1.0f / Det;
else
invDet = 0.0f;
S *= invDet;
T *= invDet;
sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0) +
T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
}
}
else
{
if(closest)
VECCOPY2D(closest, v1);
return sdist[1];
}
}
/*
* calculates the closest point between point-tri (2D)
* returns that tri must be right-handed
* Returns square distance
*/
static float closest_point_in_tri2D(const float point[2], /*const*/ float tri[3][2], float closest[2])
{
float edge_di[2];
float v_point[2];
float proj[2]; //point projected over edge-dir, edge-normal (witouth normalized edge)
const float *v0 = tri[2], *v1;
float edge_slen, d; //edge squared length
int i;
const float *nearest_vertex = NULL;
float tmp0, tmp1, numer, denom;
//for each edge
for(i=0, v0=tri[2], v1=tri[0]; i < 3; v0=tri[i++], v1=tri[i])
{
VECSUB2D(edge_di, v1, v0);
VECSUB2D(v_point, point, v0);
proj[1] = v_point[0]*edge_di[1] - v_point[1]*edge_di[0]; //dot product with edge normal
//point inside this edge
if(proj[1] < 0)
continue;
proj[0] = v_point[0]*edge_di[0] + v_point[1]*edge_di[1];
//closest to this edge is v0
if(proj[0] < 0)
if ( S < 0.0f ) // Region 2
{
if(nearest_vertex == NULL || nearest_vertex == v0)
nearest_vertex = v0;
tmp0 = A01 + B0;
tmp1 = A11 + B1;
if ( tmp1 > tmp0 )
{
numer = tmp1 - tmp0;
denom = A00 - 2.0f * A01 + A11;
if ( numer >= denom )
{
S = 1.0f;
T = 0.0f;
sqrDist = A00 + 2.0f * B0 + C;
lv = 1;
}
else
{
if(fabs(denom) > FLT_EPSILON)
S = numer / denom;
else
S = 0.0f;
T = 1.0f - S;
sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
le = 2;
}
}
else
{
//choose nearest
return choose_nearest(nearest_vertex, v0, point, closest);
S = 0.0f;
if ( tmp1 <= 0.0f )
{
T = 1.0f;
sqrDist = A11 + 2.0f * B1 + C;
lv = 2;
}
else if ( B1 >= 0.0f )
{
T = 0.0f;
sqrDist = C;
lv = 0;
}
else
{
if(fabs(A11) > FLT_EPSILON)
T = -B1 / A11;
else
T = 0.0f;
sqrDist = B1 * T + C;
le = 1;
}
}
i++; //We can skip next edge
continue;
}
edge_slen = edge_di[0]*edge_di[0] + edge_di[1]*edge_di[1]; //squared edge len
//closest to this edge is v1
if(proj[0] > edge_slen)
else if ( T < 0.0f ) // Region 6
{
if(nearest_vertex == NULL || nearest_vertex == v1)
nearest_vertex = v1;
tmp0 = A01 + B1;
tmp1 = A00 + B0;
if ( tmp1 > tmp0 )
{
numer = tmp1 - tmp0;
denom = A00 - 2.0f * A01 + A11;
if ( numer >= denom )
{
T = 1.0f;
S = 0.0f;
sqrDist = A11 + 2.0f * B1 + C;
lv = 2;
}
else
{
if(fabs(denom) > FLT_EPSILON)
T = numer / denom;
else
T = 0.0f;
S = 1.0f - T;
sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
le = 2;
}
}
else
{
return choose_nearest(nearest_vertex, v1, point, closest);
T = 0.0f;
if ( tmp1 <= 0.0f )
{
S = 1.0f;
sqrDist = A00 + 2.0f * B0 + C;
lv = 1;
}
else if ( B0 >= 0.0f )
{
S = 0.0f;
sqrDist = C;
lv = 0;
}
else
{
if(fabs(A00) > FLT_EPSILON)
S = -B0 / A00;
else
S = 0.0f;
sqrDist = B0 * S + C;
le = 0;
}
}
}
else // Region 1
{
numer = A11 + B1 - A01 - B0;
if ( numer <= 0.0f )
{
S = 0.0f;
T = 1.0f;
sqrDist = A11 + 2.0f * B1 + C;
lv = 2;
}
else
{
denom = A00 - 2.0f * A01 + A11;
if ( numer >= denom )
{
S = 1.0f;
T = 0.0f;
sqrDist = A00 + 2.0f * B0 + C;
lv = 1;
}
else
{
if(fabs(denom) > FLT_EPSILON)
S = numer / denom;
else
S = 0.0f;
T = 1.0f - S;
sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
le = 2;
}
}
continue;
}
//nearest is on this edge
d= proj[1] / edge_slen;
closest[0] = point[0] - edge_di[1] * d;
closest[1] = point[1] + edge_di[0] * d;
return proj[1]*proj[1]/edge_slen;
}
if(nearest_vertex)
// Account for numerical round-off error
if ( sqrDist < FLT_EPSILON )
sqrDist = 0.0f;
{
VECSUB2D(v_point, nearest_vertex, point);
VECCOPY2D(closest, nearest_vertex);
return v_point[0]*v_point[0] + v_point[1]*v_point[1];
float w[3], x[3], y[3], z[3];
VECCOPY(w, v0);
VECCOPY(x, e0);
VecMulf(x, S);
VECCOPY(y, e1);
VecMulf(y, T);
VECADD(z, w, x);
VECADD(z, z, y);
//VECSUB(d, p, z);
VECCOPY(nearest, z);
// d = p - ( v0 + S * e0 + T * e1 );
}
else
{
VECCOPY(closest, point); //point is already inside
return 0.0f;
}
}
*v = lv;
*e = le;
/*
* Returns the square of the minimum distance between the point and a triangle surface
* If nearest is not NULL the nearest surface point is written on it
*/
static float nearest_point_in_tri_surface(const float *point, const float *v0, const float *v1, const float *v2, float *nearest)
{
//Lets solve the 2D problem (closest point-tri)
float normal_dist, plane_sdist, plane_offset;
float du[3], dv[3], dw[3]; //orthogonal axis (du=(v0->v1), dw=plane normal)
float p_2d[2], tri_2d[3][2], nearest_2d[2];
CalcNormFloat((float*)v0, (float*)v1, (float*)v2, dw);
//point-plane distance and calculate axis
normal_dist = point_plane_distance(point, v0, dw);
// OPTIMIZATION
// if we are only interested in nearest distance if its closer than some distance already found
// we can:
// if(normal_dist*normal_dist >= best_dist_so_far) return FLOAT_MAX;
//
VECSUB(du, v1, v0);
Normalize(du);
Crossf(dv, dw, du);
plane_offset = INPR(v0, dw);
//project stuff to 2d
tri_2d[0][0] = INPR(du, v0);
tri_2d[0][1] = INPR(dv, v0);
tri_2d[1][0] = INPR(du, v1);
tri_2d[1][1] = INPR(dv, v1);
tri_2d[2][0] = INPR(du, v2);
tri_2d[2][1] = INPR(dv, v2);
p_2d[0] = INPR(du, point);
p_2d[1] = INPR(dv, point);
//we always have a right-handed tri
//this should always happen because of the way normal is calculated
plane_sdist = closest_point_in_tri2D(p_2d, tri_2d, nearest_2d);
//project back to 3d
if(nearest)
{
nearest[0] = du[0]*nearest_2d[0] + dv[0] * nearest_2d[1] + dw[0] * plane_offset;
nearest[1] = du[1]*nearest_2d[0] + dv[1] * nearest_2d[1] + dw[1] * plane_offset;
nearest[2] = du[2]*nearest_2d[0] + dv[2] * nearest_2d[1] + dw[2] * plane_offset;
}
return plane_sdist + normal_dist*normal_dist;
return sqrDist;
}
@@ -268,22 +419,15 @@ static void mesh_faces_nearest_point(void *userdata, int index, const float *co,
do
{
float nearest_tmp[3], dist;
float vec[3][3];
int vertex, edge;
// only insert valid triangles / quads with area > 0
VECSUB(vec[0], t2, t1);
VECSUB(vec[1], t0, t1);
Crossf(vec[2], vec[0], vec[1]);
if(INPR(vec[2], vec[2]) >= FLT_EPSILON)
dist = nearest_point_in_tri_surface(t0, t1, t2, co, &vertex, &edge, nearest_tmp);
if(dist < nearest->dist)
{
dist = nearest_point_in_tri_surface(co,t0, t1, t2, nearest_tmp);
if(dist < nearest->dist)
{
nearest->index = index;
nearest->dist = dist;
VECCOPY(nearest->co, nearest_tmp);
CalcNormFloat((float*)t0, (float*)t1, (float*)t2, nearest->no); //TODO.. (interpolate normals from the vertexs coordinates?
}
nearest->index = index;
nearest->dist = dist;
VECCOPY(nearest->co, nearest_tmp);
CalcNormFloat(t0, t1, t2, nearest->no);
}
t1 = t2;

View File

@@ -405,7 +405,7 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
//Now, everything is ready to project the vertexs!
//#pragma omp parallel for private(i,hit) schedule(static)
#pragma omp parallel for private(i,hit) schedule(static)
for(i = 0; i<calc->numVerts; ++i)
{
float *co = calc->vertexCos[i];

View File

@@ -69,6 +69,8 @@ variables on the UI for now
#include "BLI_blenlib.h"
#include "BLI_arithb.h"
#include "BLI_ghash.h"
#include "BLI_threads.h"
#include "BKE_curve.h"
#include "BKE_effect.h"
#include "BKE_global.h"
@@ -118,6 +120,20 @@ typedef struct SBScratch {
float aabbmin[3],aabbmax[3];
}SBScratch;
typedef struct SB_thread_context{
Object *ob;
float forcetime;
float timenow;
int ifirst;
int ilast;
ListBase *do_effector;
int do_deflector;
float fieldfactor;
float windfactor;
int nr;
int tot;
}SB_thread_context;
#define NLF_BUILD 1
#define NLF_SOLVE 2
@@ -1514,17 +1530,15 @@ int sb_detect_edge_collisionCached(float edge_v1[3],float edge_v2[3],float *damp
void scan_for_ext_spring_forces(Object *ob,float timenow)
void _scan_for_ext_spring_forces(Object *ob,float timenow,int ifirst,int ilast, struct ListBase *do_effector)
{
SoftBody *sb = ob->soft;
ListBase *do_effector;
int a;
float damp;
float feedback[3];
do_effector= pdInitEffectors(ob,NULL);
if (sb && sb->totspring){
for(a=0; a<sb->totspring; a++) {
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;
@@ -1584,9 +1598,88 @@ void scan_for_ext_spring_forces(Object *ob,float timenow)
}
}
}
if(do_effector)
pdEndEffectors(do_effector);
}
void scan_for_ext_spring_forces(Object *ob,float timenow)
{
SoftBody *sb = ob->soft;
ListBase *do_effector= NULL;
do_effector= pdInitEffectors(ob,NULL);
if (sb){
_scan_for_ext_spring_forces(ob,timenow,0,sb->totspring,do_effector);
}
if(do_effector)
pdEndEffectors(do_effector);
}
void *exec_scan_for_ext_spring_forces(void *data)
{
SB_thread_context *pctx = (SB_thread_context*)data;
_scan_for_ext_spring_forces(pctx->ob,pctx->timenow,pctx->ifirst,pctx->ilast,pctx->do_effector);
return 0;
}
void sb_sfesf_threads_run(struct Object *ob, float timenow,int totsprings,int *ptr_to_break_func())
{
ListBase *do_effector = NULL;
ListBase threads;
SB_thread_context *sb_threads;
int i, totthread,left,dec;
int lowsprings =10; /* wild guess .. may increase with better thread management 'above' or even be UI option sb->spawn_cf_threads_nopts */
do_effector= pdInitEffectors(ob,NULL);
/* figure the number of threads while preventing pretty pointless threading overhead */
if(totsprings < lowsprings) {totthread=1;}
else{
if(G.scene->r.mode & R_FIXED_THREADS)
totthread= G.scene->r.threads;
else
totthread= BLI_system_thread_count();
}
/*left to do--> what if we got zillions of CPUs running but 'totsprings' tasks to spread*/
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].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 = 0;// 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);
if(do_effector)
pdEndEffectors(do_effector);
}
/* --- the spring external section*/
int choose_winner(float*w, float* pos,float*a,float*b,float*c,float*ca,float*cb,float*cc)
@@ -2023,31 +2116,277 @@ static void sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float fo
}
static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags)
/* 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 */
int _softbody_calc_forces_slice_in_a_thread(Object *ob, float forcetime, float timenow,int ifirst,int ilast,int *ptr_to_break_func(),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);
VecSubf(def, bp->pos, obp->pos);
/* rather check the AABBoxes before ever calulating the real distance */
/* mathematically it is completly nuts, but performace is pretty much (3) times faster */
if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
distance = Normalize(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 ;
VecMidf(velcenter, bp->vec, obp->vec);
VecSubf(dvel,velcenter,bp->vec);
VecMulf(dvel,sb->nodemass);
Vec3PlusStVec(bp->force,f*(1.0f-sb->balldamp),def);
Vec3PlusStVec(bp->force,sb->balldamp,dvel);
/* exploit force(a,b) == -force(b,a) part2/2 */
VecSubf(dvel,velcenter,obp->vec);
VecMulf(dvel,sb->nodemass);
Vec3PlusStVec(obp->force,sb->balldamp,dvel);
Vec3PlusStVec(obp->force,-f*(1.0f-sb->balldamp),def);
}
}
}
}
/* naive ball self collision done */
if(bp->goal < SOFTGOALSNAP){ /* ommit 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;
VecSubf(auxvect,bp->pos,bp->origT);
ks = 1.0f/(1.0f- bp->goal*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*/
VecSubf(velgoal,bp->origS, bp->origE);
kd = sb->goalfrict * sb_fric_force_scale(ob) ;
VecAddf(auxvect,velgoal,bp->vec);
if (forcetime > 0.0 ) { /* 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 (sb){
float gravity = sb->grav * sb_grav_force_scale(ob);
bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
}
/* particle field & vortex */
if(do_effector) {
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 */
pdDoEffectors(do_effector, bp->pos, force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);
/* apply forcefield*/
VecMulf(force,fieldfactor* eval_sb_fric_force_scale);
VECADD(bp->force, 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->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;
float kd = 1.0f;
if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
VECSUB(cfforce,bp->vec,vel);
Vec3PlusStVec(bp->force,-cf*50.0f,cfforce);
Vec3PlusStVec(bp->force,kd,defforce);
}
}
/* ---cached collision targets */
/* +++springs */
iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
if(ob->softflag & OB_SB_EDGES) {
if (sb->bspring){ /* spring list exists at all ? */
int b;
BodySpring *bs;
for(b=bp->nofsprings;b>0;b--){
bs = sb->bspring + bp->springs[b-1];
if (do_springcollision || do_aero){
VecAddf(bp->force,bp->force,bs->ext_force);
if (bs->flag & BSF_INTERSECT)
bp->choke = bs->cf;
}
// sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
sb_spring_force(ob,ilast-bb,bs,iks,forcetime,0);
}/* loop springs */
}/* existing spring list */
}/*any edges*/
/* ---springs */
}/*omit on snap */
}/*loop all bp's*/
return 0; /*done fine*/
}
void *exec_softbody_calc_forces(void *data)
{
SB_thread_context *pctx = (SB_thread_context*)data;
_softbody_calc_forces_slice_in_a_thread(pctx->ob,pctx->forcetime,pctx->timenow,pctx->ifirst,pctx->ilast,NULL,pctx->do_effector,pctx->do_deflector,pctx->fieldfactor,pctx->windfactor);
return 0;
}
void sb_cf_threads_run(struct Object *ob, float forcetime, float timenow,int totpoint,int *ptr_to_break_func(),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 =10; /* 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 */
if(totpoint < lowpoints) {totthread=1;}
else{
if(G.scene->r.mode & R_FIXED_THREADS)
totthread= G.scene->r.threads;
else
totthread= BLI_system_thread_count();
}
/*left to do--> what if we got zillions of CPUs running but 'totpoint' tasks to spread*/
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].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(Object *ob, float forcetime, float timenow, int nl_flags)
{
/* 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;
BodySpring *bs;
ListBase *do_effector;
float iks, ks, kd, gravity;
float iks, gravity;
float fieldfactor = 1000.0f, windfactor = 250.0f;
float tune = sb->ballstiff;
int a, b, do_deflector,do_selfcollision,do_springcollision,do_aero;
/* jacobian
NLboolean success;
if(nl_flags){
nlBegin(NL_SYSTEM);
nlBegin(NL_MATRIX);
}
*/
int do_deflector,do_selfcollision,do_springcollision,do_aero;
gravity = sb->grav * sb_grav_force_scale(ob);
@@ -2060,7 +2399,9 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
iks = 1.0f/(1.0f-sb->inspring)-1.0f ;/* inner spring constants function */
bproot= sb->bpoint; /* need this for proper spring addressing */
if (do_springcollision || do_aero) scan_for_ext_spring_forces(ob,timenow);
if (do_springcollision || do_aero)
sb_sfesf_threads_run(ob,timenow,sb->totspring,NULL);
/* after spring scan because it uses Effoctors too */
do_effector= pdInitEffectors(ob,NULL);
@@ -2069,34 +2410,102 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
}
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
/* clear forces accumulator */
bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
//int op =3*sb->totpoint;
/* dF/dV = v */
/* jacobioan
nlMatrixAdd(op+ia,ia,-forcetime);
nlMatrixAdd(op+ia+1,ia+1,-forcetime);
nlMatrixAdd(op+ia+2,ia+2,-forcetime);
nlMatrixAdd(ia,ia,1);
nlMatrixAdd(ia+1,ia+1,1);
nlMatrixAdd(ia+2,ia+2,1);
sb_cf_threads_run(ob,forcetime,timenow,sb->totpoint,NULL,do_effector,do_deflector,fieldfactor,windfactor);
nlMatrixAdd(op+ia,op+ia,1);
nlMatrixAdd(op+ia+1,op+ia+1,1);
nlMatrixAdd(op+ia+2,op+ia+2,1);
*/
/* finally add forces caused by face collision */
if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
/* finish matrix and solve */
if(do_effector) pdEndEffectors(do_effector);
}
static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int nl_flags)
{
/* redirection to the new threaded Version */
if (G.rt !=16){
softbody_calc_forcesEx(ob, forcetime, timenow, nl_flags);
return;
}
else{
/* so the following will die */
/* |||||||||||||||||||||||||| */
/* VVVVVVVVVVVVVVVVVVVVVVVVVV */
/* 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;
BodySpring *bs;
ListBase *do_effector;
float iks, ks, kd, gravity;
float fieldfactor = 1000.0f, windfactor = 250.0f;
float tune = sb->ballstiff;
int a, b, do_deflector,do_selfcollision,do_springcollision,do_aero;
/* jacobian
NLboolean success;
if(nl_flags){
nlBegin(NL_SYSTEM);
nlBegin(NL_MATRIX);
}
*/
gravity = sb->grav * sb_grav_force_scale(ob);
/* check conditions for various options */
do_deflector= query_external_colliders(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 */
if (do_springcollision || do_aero) scan_for_ext_spring_forces(ob,timenow);
/* after spring scan because it uses Effoctors too */
do_effector= pdInitEffectors(ob,NULL);
if (do_deflector) {
float defforce[3];
do_deflector = sb_detect_aabb_collisionCached(defforce,ob->lay,ob,timenow);
}
/* naive ball self collision */
/* needs to be done if goal snaps or not */
if(do_selfcollision){
int attached;
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
/* clear forces accumulator */
bp->force[0]= bp->force[1]= bp->force[2]= 0.0;
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
//int op =3*sb->totpoint;
/* dF/dV = v */
/* jacobioan
nlMatrixAdd(op+ia,ia,-forcetime);
nlMatrixAdd(op+ia+1,ia+1,-forcetime);
nlMatrixAdd(op+ia+2,ia+2,-forcetime);
nlMatrixAdd(ia,ia,1);
nlMatrixAdd(ia+1,ia+1,1);
nlMatrixAdd(ia+2,ia+2,1);
nlMatrixAdd(op+ia,op+ia,1);
nlMatrixAdd(op+ia+1,op+ia+1,1);
nlMatrixAdd(op+ia+2,op+ia+2,1);
*/
}
/* naive ball self collision */
/* needs to be done if goal snaps or not */
if(do_selfcollision){
int attached;
BodyPoint *obp;
int c,b;
float velcenter[3],dvel[3],def[3];
@@ -2104,7 +2513,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
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);
@@ -2114,7 +2523,7 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
/* mathematically it is completly nuts, but performace is pretty much (3) times faster */
if ((ABS(def[0]) > compare) || (ABS(def[1]) > compare) || (ABS(def[2]) > compare)) continue;
distance = Normalize(def);
distance = Normalize(def);
if (distance < compare ){
/* exclude body points attached with a spring */
attached = 0;
@@ -2146,8 +2555,8 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
/* exploit force(a,b) == -force(b,a) part1/2 */
//dfdx_goal(ic,ic,op,mpos);
//dfdv_goal(ic,ic,mvel);
/*TODO sit down an X-out the true jacobian entries*/
/*well does not make to much sense because the eigenvalues
of the jacobian go negative; and negative eigenvalues
@@ -2169,236 +2578,239 @@ static void softbody_calc_forces(Object *ob, float forcetime, float timenow, int
}
}
}
}
/* naive ball self collision done */
}
/* naive ball self collision done */
if(bp->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */
float auxvect[3];
float velgoal[3];
if(bp->goal < SOFTGOALSNAP){ /* ommit this bp when it snaps */
float auxvect[3];
float velgoal[3];
/* do goal stuff */
if(ob->softflag & OB_SB_GOAL) {
/* true elastic goal */
VecSubf(auxvect,bp->pos,bp->origT);
ks = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
bp->force[0]+= -ks*(auxvect[0]);
bp->force[1]+= -ks*(auxvect[1]);
bp->force[2]+= -ks*(auxvect[2]);
/* do goal stuff */
if(ob->softflag & OB_SB_GOAL) {
/* true elastic goal */
VecSubf(auxvect,bp->pos,bp->origT);
ks = 1.0f/(1.0f- bp->goal*sb->goalspring)-1.0f ;
bp->force[0]+= -ks*(auxvect[0]);
bp->force[1]+= -ks*(auxvect[1]);
bp->force[2]+= -ks*(auxvect[2]);
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
//int op =3*(sb->totpoint);
/* depending on my pos */
//dfdx_goal(ia,ia,op,ks*forcetime);
}
/* calulate damping forces generated by goals*/
VecSubf(velgoal,bp->origS, bp->origE);
kd = sb->goalfrict * sb_fric_force_scale(ob) ;
VecAddf(auxvect,velgoal,bp->vec);
if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */
bp->force[0]-= kd * (auxvect[0]);
bp->force[1]-= kd * (auxvect[1]);
bp->force[2]-= kd * (auxvect[2]);
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
Normalize(auxvect);
/* depending on my vel */
//dfdv_goal(ia,ia,kd*forcetime);
//int op =3*(sb->totpoint);
/* depending on my pos */
//dfdx_goal(ia,ia,op,ks*forcetime);
}
/* calulate damping forces generated by goals*/
VecSubf(velgoal,bp->origS, bp->origE);
kd = sb->goalfrict * sb_fric_force_scale(ob) ;
VecAddf(auxvect,velgoal,bp->vec);
if (forcetime > 0.0 ) { /* make sure friction does not become rocket motor on time reversal */
bp->force[0]-= kd * (auxvect[0]);
bp->force[1]-= kd * (auxvect[1]);
bp->force[2]-= kd * (auxvect[2]);
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
Normalize(auxvect);
/* depending on my vel */
//dfdv_goal(ia,ia,kd*forcetime);
}
}
else {
bp->force[0]-= kd * (velgoal[0] - bp->vec[0]);
bp->force[1]-= kd * (velgoal[1] - bp->vec[1]);
bp->force[2]-= kd * (velgoal[2] - bp->vec[2]);
}
}
/* done goal stuff */
/* gravitation */
bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
//bp->force[1]-= gravity*sb->nodemass; /* individual mass of node here */
/* particle field & vortex */
if(do_effector) {
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 */
pdDoEffectors(do_effector, bp->pos, force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);
/* apply forcefield*/
VecMulf(force,fieldfactor* eval_sb_fric_force_scale);
VECADD(bp->force, 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->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 */
bp->force[2]-= gravity*sb->nodemass; /* individual mass of node here */
//bp->force[1]-= gravity*sb->nodemass; /* individual mass of node here */
/* BP friction in media (not) moving*/
kd= sb->mediafrict* sb_fric_force_scale(ob);
/* assume it to be proportional to actual velocity */
bp->force[0]-= bp->vec[0]*kd;
bp->force[1]-= bp->vec[1]*kd;
bp->force[2]-= bp->vec[2]*kd;
/* friction in media done */
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
/* da/dv = */
/* particle field & vortex */
if(do_effector) {
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 */
pdDoEffectors(do_effector, bp->pos, force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);
/* apply forcefield*/
VecMulf(force,fieldfactor* eval_sb_fric_force_scale);
VECADD(bp->force, bp->force, force);
/* BP friction in moving media */
kd= sb->mediafrict* eval_sb_fric_force_scale;
bp->force[0] -= kd * (bp->vec[0] + windfactor*speed[0]/eval_sb_fric_force_scale);
bp->force[1] -= kd * (bp->vec[1] + windfactor*speed[1]/eval_sb_fric_force_scale);
bp->force[2] -= kd * (bp->vec[2] + windfactor*speed[2]/eval_sb_fric_force_scale);
/* now we'll have nice centrifugal effect for vortex */
}
else {
/* BP friction in media (not) moving*/
kd= sb->mediafrict* sb_fric_force_scale(ob);
/* assume it to be proportional to actual velocity */
bp->force[0]-= bp->vec[0]*kd;
bp->force[1]-= bp->vec[1]*kd;
bp->force[2]-= bp->vec[2]*kd;
/* friction in media done */
if(nl_flags & NLF_BUILD){
//int ia =3*(sb->totpoint-a);
/* da/dv = */
// nlMatrixAdd(ia,ia,forcetime*kd);
// nlMatrixAdd(ia+1,ia+1,forcetime*kd);
// nlMatrixAdd(ia+2,ia+2,forcetime*kd);
}
}
/* +++cached collision targets */
bp->choke = 0.0f;
bp->choke2 = 0.0f;
bp->flag &= ~SBF_DOFUZZY;
if(do_deflector) {
float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion;
kd = 1.0f;
if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
if ((!nl_flags)&&(intrusion < 0.0f)){
/*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 beeing out of the collider
in heun step No1 and leave the heun step No2 adapt to it
so we kind of introduced a implicit solver for this case
*/
Vec3PlusStVec(bp->pos,-intrusion,facenormal);
sb->scratch->flag |= SBF_DOFUZZY;
bp->flag |= SBF_DOFUZZY;
bp->choke = sb->choke*0.01f;
// nlMatrixAdd(ia,ia,forcetime*kd);
// nlMatrixAdd(ia+1,ia+1,forcetime*kd);
// nlMatrixAdd(ia+2,ia+2,forcetime*kd);
}
else{
}
/* +++cached collision targets */
bp->choke = 0.0f;
bp->choke2 = 0.0f;
bp->flag &= ~SBF_DOFUZZY;
if(do_deflector) {
float cfforce[3],defforce[3] ={0.0f,0.0f,0.0f}, vel[3] = {0.0f,0.0f,0.0f}, facenormal[3], cf = 1.0f,intrusion;
kd = 1.0f;
if (sb_deflect_face(ob,bp->pos,facenormal,defforce,&cf,timenow,vel,&intrusion)){
if ((!nl_flags)&&(intrusion < 0.0f)){
/*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 beeing out of the collider
in heun step No1 and leave the heun step No2 adapt to it
so we kind of introduced a implicit solver for this case
*/
Vec3PlusStVec(bp->pos,-intrusion,facenormal);
sb->scratch->flag |= SBF_DOFUZZY;
bp->flag |= SBF_DOFUZZY;
bp->choke = sb->choke*0.01f;
}
else{
VECSUB(cfforce,bp->vec,vel);
Vec3PlusStVec(bp->force,-cf*50.0f,cfforce);
}
Vec3PlusStVec(bp->force,kd,defforce);
if (nl_flags & NLF_BUILD){
// int ia =3*(sb->totpoint-a);
// int op =3*sb->totpoint;
//dfdx_goal(ia,ia,op,mpos); // don't do unless you know
//dfdv_goal(ia,ia,-cf);
}
}
}
/* ---cached collision targets */
/* +++springs */
if(ob->softflag & OB_SB_EDGES) {
if (sb->bspring){ /* spring list exists at all ? */
for(b=bp->nofsprings;b>0;b--){
bs = sb->bspring + bp->springs[b-1];
if (do_springcollision || do_aero){
VecAddf(bp->force,bp->force,bs->ext_force);
if (bs->flag & BSF_INTERSECT)
bp->choke = bs->cf;
}
Vec3PlusStVec(bp->force,kd,defforce);
if (nl_flags & NLF_BUILD){
// int ia =3*(sb->totpoint-a);
// int op =3*sb->totpoint;
//dfdx_goal(ia,ia,op,mpos); // don't do unless you know
//dfdv_goal(ia,ia,-cf);
}
// sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,nl_flags);
}/* loop springs */
}/* existing spring list */
}/*any edges*/
/* ---springs */
}/*omit on snap */
}/*loop all bp's*/
}
}
/* ---cached collision targets */
/* +++springs */
if(ob->softflag & OB_SB_EDGES) {
if (sb->bspring){ /* spring list exists at all ? */
for(b=bp->nofsprings;b>0;b--){
bs = sb->bspring + bp->springs[b-1];
if (do_springcollision || do_aero){
VecAddf(bp->force,bp->force,bs->ext_force);
if (bs->flag & BSF_INTERSECT)
bp->choke = bs->cf;
}
// sb_spring_force(Object *ob,int bpi,BodySpring *bs,float iks,float forcetime,int nl_flags)
// rather remove nl_falgs from code .. will make things a lot cleaner
sb_spring_force(ob,sb->totpoint-a,bs,iks,forcetime,0);
}/* loop springs */
}/* existing spring list */
}/*any edges*/
/* ---springs */
}/*omit on snap */
}/*loop all bp's*/
/* finally add forces caused by face collision */
if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
/* finish matrix and solve */
/* finally add forces caused by face collision */
if (ob->softflag & OB_SB_FACECOLL) scan_for_ext_face_forces(ob,timenow);
/* finish matrix and solve */
#if (0) // remove onl linking for now .. still i am not sure .. the jacobian can be usefull .. so keep that BM
if(nl_flags & NLF_SOLVE){
//double sct,sst=PIL_check_seconds_timer();
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
int iv =3*(sb->totpoint-a);
int ip =3*(2*sb->totpoint-a);
int n;
for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);}
for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);}
}
nlEnd(NL_MATRIX);
nlEnd(NL_SYSTEM);
if ((G.rt >0) && (nl_flags & NLF_BUILD))
{
printf("####MEE#####\n");
nlPrintMatrix();
}
success= nlSolveAdvanced(NULL, 1);
// nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */
if(success){
float f;
int index =0;
/* for debug purpose .. anyhow cropping B vector looks like working */
if (G.rt >0)
for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
printf("(%f ",f);index++;
f=nlGetVariable(0,index);
printf("%f ",f);index++;
f=nlGetVariable(0,index);
printf("%f)",f);index++;
}
index =0;
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
bp->impdv[0] = f; index++;
f=nlGetVariable(0,index);
bp->impdv[1] = f; index++;
f=nlGetVariable(0,index);
bp->impdv[2] = f; index++;
}
/*
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
bp->impdx[0] = f; index++;
f=nlGetVariable(0,index);
bp->impdx[1] = f; index++;
f=nlGetVariable(0,index);
bp->impdx[2] = f; index++;
}
*/
}
else{
printf("Matrix inversion failed \n");
if(nl_flags & NLF_SOLVE){
//double sct,sst=PIL_check_seconds_timer();
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
VECCOPY(bp->impdv,bp->force);
int iv =3*(sb->totpoint-a);
int ip =3*(2*sb->totpoint-a);
int n;
for (n=0;n<3;n++) {nlRightHandSideSet(0, iv+n, bp->force[0+n]);}
for (n=0;n<3;n++) {nlRightHandSideSet(0, ip+n, bp->vec[0+n]);}
}
nlEnd(NL_MATRIX);
nlEnd(NL_SYSTEM);
if ((G.rt == 32) && (nl_flags & NLF_BUILD))
{
printf("####MEE#####\n");
nlPrintMatrix();
}
}
success= nlSolveAdvanced(NULL, 1);
//sct=PIL_check_seconds_timer();
//if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name);
}
/* cleanup */
// nlPrintMatrix(); /* for debug purpose .. anyhow cropping B vector looks like working */
if(success){
float f;
int index =0;
/* for debug purpose .. anyhow cropping B vector looks like working */
if (G.rt ==32)
for(a=2*sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
printf("(%f ",f);index++;
f=nlGetVariable(0,index);
printf("%f ",f);index++;
f=nlGetVariable(0,index);
printf("%f)",f);index++;
}
index =0;
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
bp->impdv[0] = f; index++;
f=nlGetVariable(0,index);
bp->impdv[1] = f; index++;
f=nlGetVariable(0,index);
bp->impdv[2] = f; index++;
}
/*
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
f=nlGetVariable(0,index);
bp->impdx[0] = f; index++;
f=nlGetVariable(0,index);
bp->impdx[1] = f; index++;
f=nlGetVariable(0,index);
bp->impdx[2] = f; index++;
}
*/
}
else{
printf("Matrix inversion failed \n");
for(a=sb->totpoint, bp= sb->bpoint; a>0; a--, bp++) {
VECCOPY(bp->impdv,bp->force);
}
}
//sct=PIL_check_seconds_timer();
//if (sct-sst > 0.01f) printf(" implicit solver time %f %s \r",sct-sst,ob->id.name);
}
/* cleanup */
#endif
if(do_effector) pdEndEffectors(do_effector);
if(do_effector) pdEndEffectors(do_effector);
}
}
static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *err, int mid_flags)
{
/* time evolution */
@@ -2458,7 +2870,7 @@ static void softbody_apply_forces(Object *ob, float forcetime, int mode, float *
/* x(t + dt) = x(t) + v(t~) * dt */
VecMulf(dx,forcetime);
/* the freezer */
/* the freezer coming sooner or later */
/*
if ((Inpf(dx,dx)<freezeloc )&&(Inpf(bp->force,bp->force)<freezeforce )){
bp->frozen /=2;
@@ -3529,6 +3941,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
* 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 */
if (sb->minloops > 0) forcetimemax = 1.0f / sb->minloops;
@@ -3546,13 +3959,13 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
sb->scratch->flag &= ~SBF_DOFUZZY;
/* do predictive euler step */
softbody_calc_forces(ob, forcetime,timedone/dtime,0);
softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags);
softbody_apply_forces(ob, forcetime, 1, NULL,mid_flags);
/* crop new slope values to do averaged slope step */
softbody_calc_forces(ob, forcetime,timedone/dtime,0);
softbody_apply_forces(ob, forcetime, 2, &err,mid_flags);
softbody_apply_forces(ob, forcetime, 2, &err,mid_flags);
softbody_apply_goalsnap(ob);
if (err > SoftHeunTol) { /* error needs to be scaled to some quantity */
@@ -3603,7 +4016,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
// if(G.f & G_DEBUG){
if(sb->solverflags & SBSO_MONITOR ){
if (loops > HEUNWARNLIMIT) /* monitor high loop counts */
printf("\r needed %d steps/frame ",loops);
printf("\r needed %d steps/frame",loops);
}
}
@@ -3627,7 +4040,7 @@ static void softbody_step(Object *ob, SoftBody *sb, float dtime)
if(sb->solverflags & SBSO_MONITOR ){
sct=PIL_check_seconds_timer();
if (sct-sst > 0.5f) printf(" solver time %f %s \r",sct-sst,ob->id.name);
if (sct-sst > 0.5f) printf(" solver time %f sec %s \n",sct-sst,ob->id.name);
}
}

View File

@@ -105,6 +105,10 @@
#define EXPP_TEX_LACUNARITY_MAX 6.0f
#define EXPP_TEX_OCTS_MIN 0.0f
#define EXPP_TEX_OCTS_MAX 8.0f
#define EXPP_TEX_OFST_MIN 0.0f
#define EXPP_TEX_OFST_MAX 6.0f
#define EXPP_TEX_GAIN_MIN 0.0f
#define EXPP_TEX_GAIN_MAX 6.0f
#define EXPP_TEX_ISCALE_MIN 0.0f
#define EXPP_TEX_ISCALE_MAX 10.0f
#define EXPP_TEX_EXP_MIN 0.010f
@@ -430,6 +434,8 @@ GETFUNC( getNoiseDepth );
GETFUNC( getNoiseSize );
GETFUNC( getNoiseType );
GETFUNC( getOcts );
GETFUNC( getOffset );
GETFUNC( getGain );
GETFUNC( getRepeat );
GETFUNC( getRGBCol );
GETFUNC( getSType );
@@ -478,6 +484,8 @@ SETFUNC( setNoiseDepth );
SETFUNC( setNoiseSize );
SETFUNC( setNoiseType );
SETFUNC( setOcts );
SETFUNC( setOffset );
SETFUNC( setGain );
SETFUNC( setRepeat );
SETFUNC( setRGBCol );
SETFUNC( setSType );
@@ -646,6 +654,14 @@ static PyGetSetDef BPy_Texture_getseters[] = {
(getter)Texture_getLacunarity, (setter)Texture_setLacunarity,
"Gap between succesive frequencies (for Musgrave textures)",
NULL},
{"offset",
(getter)Texture_getOffset, (setter)Texture_setOffset,
"Fractal offset (for Musgrave textures)",
NULL},
{"gain",
(getter)Texture_getGain, (setter)Texture_setGain,
"Gain multiplier (for Musgrave textures)",
NULL},
{"noiseBasis",
(getter)Texture_getNoiseBasis, (setter)Texture_setNoiseBasis,
"Noise basis type (wood, stucci, marble, clouds, Musgrave, distorted noise)",
@@ -1837,6 +1853,20 @@ static int Texture_setLacunarity( BPy_Texture * self, PyObject * value )
EXPP_TEX_LACUNARITY_MAX );
}
static int Texture_setOffset( BPy_Texture * self, PyObject * value )
{
return EXPP_setFloatClamped ( value, &self->texture->mg_offset,
EXPP_TEX_OFST_MIN,
EXPP_TEX_OFST_MAX );
}
static int Texture_setGain( BPy_Texture * self, PyObject * value )
{
return EXPP_setFloatClamped ( value, &self->texture->mg_gain,
EXPP_TEX_GAIN_MIN,
EXPP_TEX_GAIN_MAX );
}
static int Texture_setOcts( BPy_Texture * self, PyObject * value )
{
return EXPP_setFloatClamped ( value, &self->texture->mg_octaves,
@@ -2168,6 +2198,16 @@ static PyObject *Texture_getOcts( BPy_Texture *self )
return PyFloat_FromDouble( self->texture->mg_octaves );
}
static PyObject *Texture_getOffset( BPy_Texture *self )
{
return PyFloat_FromDouble( self->texture->mg_offset );
}
static PyObject *Texture_getGain( BPy_Texture *self )
{
return PyFloat_FromDouble( self->texture->mg_gain );
}
static PyObject *Texture_getRepeat( BPy_Texture *self )
{
return Py_BuildValue( "(i,i)", self->texture->xrepeat,

View File

@@ -344,6 +344,12 @@ class Texture:
@ivar octs: Number of frequencies (for Musgrave textures).
Value is clamped to the range [0.0,8.0].
@type octs: float
@ivar offset: Fractal offset (for hetero terrain and multifractal Musgrave textures).
Value is clamped to the range [0.0,6.0].
@type offset: float
@ivar gain: Gain multiplier (for multifractal Musgrave textures).
Value is clamped to the range [0.0,6.0].
@type gain: float
@ivar repeat: Repetition multiplier (for image textures).
@type repeat: tuple of 2 ints
@ivar rgbCol: RGB color tuple.

View File

@@ -6720,3 +6720,4 @@ void editing_panels()
}
uiClearButLock();
}

View File

@@ -700,9 +700,7 @@ CValue* CValue::ConvertPythonToValue(PyObject* pyobj)
CValue* vallie = NULL;
PyTypeObject* type = pyobj->ob_type;
if (type == &PyList_Type)
if (PyList_Check(pyobj))
{
CListValue* listval = new CListValue();
bool error = false;
@@ -732,26 +730,25 @@ CValue* CValue::ConvertPythonToValue(PyObject* pyobj)
}
} else
if (type == &PyFloat_Type)
if (PyFloat_Check(pyobj))
{
float fl;
PyArg_Parse(pyobj,"f",&fl);
vallie = new CFloatValue(fl);
vallie = new CFloatValue( (float)PyFloat_AsDouble(pyobj) );
} else
if (type==&PyInt_Type)
if (PyInt_Check(pyobj))
{
int innie;
PyArg_Parse(pyobj,"i",&innie);
vallie = new CIntValue(innie);
vallie = new CIntValue( (int)PyInt_AS_LONG(pyobj) );
} else
if (type==&PyString_Type)
if (PyString_Check(pyobj))
{
vallie = new CStringValue(PyString_AsString(pyobj),"");
} else
if (type==&CValue::Type || type==&CListValue::Type)
if (pyobj->ob_type==&CValue::Type || pyobj->ob_type==&CListValue::Type)
{
vallie = ((CValue*) pyobj)->AddRef();
} else
{
/* return an error value from the caller */
PyErr_SetString(PyExc_TypeError, "This python value could not be assigned to a game engine property");
}
return vallie;
@@ -778,6 +775,9 @@ int CValue::_setattr(const STR_String& attr,PyObject* pyobj)
SetProperty(attr,vallie);
}
vallie->Release();
} else
{
return 1; /* ConvertPythonToValue sets the error message */
}
//PyObjectPlus::_setattr(attr,value);

View File

@@ -273,36 +273,16 @@ void SCA_PythonController::Trigger(SCA_LogicManager* logicmgr)
* break it by hand, then DECREF (which in this case
* should always ensure excdict is cleared).
*/
/* PyObject *excdict= myPyDict_Copy(m_pythondictionary);
struct _object* resultobj = PyEval_EvalCode((PyCodeObject*)m_bytecode,
excdict,
excdict
);
PyDict_Clear(excdict);
Py_DECREF(excdict);*/
#if 1
PyObject *excdict= PyDict_Copy(m_pythondictionary);
PyObject* resultobj = PyEval_EvalCode((PyCodeObject*)m_bytecode,
excdict,
excdict
);
PyDict_Clear(excdict);
Py_DECREF(excdict);
#else
PyObject* resultobj = PyEval_EvalCode((PyCodeObject*)m_bytecode,
m_pythondictionary,
m_pythondictionary
);
#endif
excdict, excdict);
if (resultobj)
{
Py_DECREF(resultobj);
} else
}
else
{
// something is wrong, tell the user what went wrong
printf("PYTHON SCRIPT ERROR:\n");
@@ -310,6 +290,11 @@ void SCA_PythonController::Trigger(SCA_LogicManager* logicmgr)
//PyRun_SimpleString(m_scriptText.Ptr());
}
// clear after PyErrPrint - seems it can be using
// something in this dictionary and crash?
PyDict_Clear(excdict);
Py_DECREF(excdict);
m_sCurrentController = NULL;
}