[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [28928] branches/hairsim/source/blender: durian hair: commit of working copy; note: will not compile easily yet, need to still add cppad to intern
Joseph Eagar
joeedh at gmail.com
Sun May 23 10:00:45 CEST 2010
Revision: 28928
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=28928
Author: joeedh
Date: 2010-05-23 10:00:44 +0200 (Sun, 23 May 2010)
Log Message:
-----------
durian hair: commit of working copy; note: will not compile easily yet, need to still add cppad to intern
Modified Paths:
--------------
branches/hairsim/source/blender/blenkernel/intern/collision.c
branches/hairsim/source/blender/blenkernel/intern/implicit.c
branches/hairsim/source/blender/blenkernel/intern/particle.c
branches/hairsim/source/blender/editors/physics/particle_edit.c
branches/hairsim/source/blender/editors/space_view3d/drawobject.c
Added Paths:
-----------
branches/hairsim/source/blender/blenkernel/intern/derivatives.cpp
Modified: branches/hairsim/source/blender/blenkernel/intern/collision.c
===================================================================
--- branches/hairsim/source/blender/blenkernel/intern/collision.c 2010-05-23 07:53:09 UTC (rev 28927)
+++ branches/hairsim/source/blender/blenkernel/intern/collision.c 2010-05-23 08:00:44 UTC (rev 28928)
@@ -66,6 +66,53 @@
#define SIGN(n) ((n) < -DBL_EPSILON)
+/*n=number, g=initial guess - (n/g + g)*0.5*/
+#define GUESS(s, g)\
+g = s*0.5;\
+ if (g > -5.0 && g < 5.0) {\
+ double g2 = 1.0 - g*0.2+DBL_EPSILON;\
+ g = guesstab[(int)ceil(g2)];\
+ } else g = g*0.5;
+
+#define FAST_SQRT(n, g, t) (t=(n/g + g)*0.5, t=(n/t + t)*0.5, t=(n/t + t)*0.5, t=(n/t + t)*0.5, t)
+/*#define normalize_fast_d(v) {\
+double _d = v[0]*v[0] + v[1]*v[1] + v[2]*v[2], _g, _t; \
+ GUESS(_d, _g);\
+ _d = FAST_SQRT(_d, _g, _t); \
+ _d = 1.0/_d; \
+ v[0]*=_d; v[1]*=_d; v[2]*=_d;\
+ }
+
+#undef normalize_fast_d*/
+#define normalize_fast_d normalize_v3_d
+
+#define GUESSTAB_LEN 8000
+static double *guesstab = NULL;
+
+static double acos_fast(double angle) {
+ return saacos_d(angle);
+
+ if (angle < -1.0) return -M_PI;
+ if (angle > 1.0) return M_PI;
+
+ return (1.0 - angle)*0.5*M_PI;
+}
+
+//#define saacos_d acos_fast
+
+static void init_guess_tab(void) {
+ int i;
+
+ if (!guesstab) {
+ guesstab = malloc(sizeof(double)*GUESSTAB_LEN);
+
+ for (i=0; i<GUESSTAB_LEN; i++) {
+ double f = (double)i / (double)GUESSTAB_LEN;
+ guesstab[i] = sqrt(f);
+ }
+ }
+}
+
static int tri_hash_comp(void *a, void *b);
/*"point-based dynamics" by mueller*/
@@ -85,7 +132,6 @@
Constraint head;
ClothVertex *v1, *v2, *v3, *p;
ClothTri *colltri;
- cdouble_t mat[4][4];
int flip;
cdouble_t mindis, no[3];
cdouble_t friction;
@@ -108,25 +154,24 @@
VECSUB(vec, v1->tx, v2->tx);
dis = len_v3_d(vec);
+
if (fabs(dis-scon->restlen) <= DBL_EPSILON) {
scon->head.result = 0.0;
return;
}
- if (scon->ignore_stretch && dis > scon->restlen) {
+ if (scon->ignore_stretch > 0 && dis > scon->restlen) {
scon->head.result = 0.0;
return;
+ } else if (scon->ignore_stretch < 0 && dis < scon->restlen) {
+ scon->head.result = 0.0;
+ return;
}
scon->head.result = (dis - scon->restlen);
mul_v3_d(vec, (1.0/dis));
- /*the idea is when you multiply grad with
- result, you get a vector straight to
- the position that satisfies the constraint.
- in this particular case, normalized vec
- works*/
copy_v3_v3_d(scon->head.grad[0], vec);
mul_v3_d(scon->head.grad[0], 1.0);
@@ -155,7 +200,7 @@
}
void solve_bending_gradient(double ip1[3], double ip2[3], double ip3[3], double grad1[3],
- double grad2[3], double grad3[3], double restangle);
+ double grad2[3], double grad3[3], double restangle);
static void eval_point_constraint(void *vcon) {
@@ -163,6 +208,7 @@
ClothVertex *v = con->head.verts[0];
cdouble_t uv[3], uv2[3], vel[3], no[3], co2[3], off[3], v1[3], v2[3], v3[3];
cdouble_t lambda;
+ cdouble_t g1[3], g2[3], g3[3];
int v1i, v2i, v3i;
VECCOPY(co2, con->head.verts[0]->tx);
@@ -178,6 +224,7 @@
VECCOPY(v3, con->colltri->verts[2]->tx);
}
normal_tri_v3_d(no, v1, v2, v3);
+ copy_v3_v3_d(con->no, no);
VECCOPY(off, no);
mul_v3_d(off, con->mindis);
@@ -186,28 +233,10 @@
add_v3_v3_d(v2, off);
add_v3_v3_d(v3, off);
- /*we don't restrict to the triangle's bounds here, the whole plane is fine*/
-#if 1
- if (isect_ray_tri_plane_v3_d(co2, no, v1, v2, v3, &lambda, uv)) {
- double cdouble_t, g1[3], g2[3], g3[3], lambda2;
-
- VECSUB(vel, v->txold, co2);
- normalize_v3_d(vel);
-
- /*probably a highly inaccurate way to approximate friction, but oh well*/
- if (isect_ray_tri_plane_v3_d(co2, vel, v1, v2, v3, &lambda2, uv2) && lambda2 < lambda*1.5) {
- VECSUB2(vel, no);
- VECMULS(vel, con->friction);
- VECADD2(no, vel);
-
- lambda += (lambda2 - lambda)*con->friction;
- }
-#endif
- //cdouble_t g1[3], g2[3], g3[3];
-
+ //if (isect_ray_tri_v3_d(co2, no, v1, v2, v3, &lambda, uv)) {
sub_v3_v3v3_d(co2, v->tx, v1);
- normalize_v3_d(co2);
- normalize_v3_d(no);
+ normalize_fast_d(co2);
+ normalize_fast_d(no);
lambda = dot_v3v3_d(co2, no);
@@ -220,12 +249,11 @@
solve_bending_gradient(v->tx, v1, co2, g1, g2, g3, M_PI/2);
VECCOPY(con->head.grad[0], g1);
- mul_v3_d(con->head.grad[0], 0.1);
- con->head.result = fabs(saacos_d(lambda));
+ con->head.result = (saacos_d(lambda) - M_PI/2);
+ //}
+ return;
//}
- return;
- }
con->head.result = 0.0;
}
@@ -265,7 +293,7 @@
VECCOPY(oco, point->txold);
VECSUB(off, oco, v1->txold);
- normalize_v3_d(off);
+ normalize_fast_d(off);
if (dot_v3v3_d(off, no) < -DBL_EPSILON) {
//con->flip = 1;
}
@@ -305,7 +333,7 @@
sub_v3_v3v3_d(vec1, v2->txconst, v1->txconst);
sub_v3_v3v3_d(vec2, v2->tx, v1->tx);
- l2 = len_v3_d(vec2);
+ l2 = len_v3_d(vec1);
normalize_v3_d(vec1);
normalize_v3_d(vec2);
@@ -321,21 +349,21 @@
sub_v3_v3v3_d(vec3, self->parent->head.verts[0]->txconst, self->parent->v1->txconst);
sub_v3_v3v3_d(vec4, self->parent->head.verts[0]->tx, self->parent->v1->tx);
- normalize_v3_d(vec3);
- normalize_v3_d(vec4);
+ normalize_fast_d(vec3);
+ normalize_fast_d(vec4);
if (self->parent->parent && !self->parent->parent->coangular) {
mul_qt_v3_d(self->parent->parent->quat, vec3);
}
- normalize_v3_d(vec3);
+ normalize_fast_d(vec3);
t = dot_v3v3_d(vec3, vec4);
- if (t < 1.0-DBL_EPSILON) {
+ if (t < 1.0-DBL_EPSILON*2) {
cross_v3_v3v3_d(self->parent->nor, vec3, vec4);
- normalize_v3_d(self->parent->nor);
+ normalize_fast_d(self->parent->nor);
- t = saacos_d(t)*0.5;
+ t = saacos_d(t)*0.15;
axis_angle_to_quat_d(self->parent->quat, self->parent->nor, t);
if (self->parent->parent) {
@@ -350,7 +378,7 @@
unit_qt_d(self->parent->quat);
self->parent->coangular = 1;
}
-
+
mul_qt_v3_d(self->parent->quat, vec1);
} else {
unit_qt_d(self->quat);
@@ -359,29 +387,28 @@
dot = dot_v3v3_d(vec1, vec2);
- /*see if we are colinear, if so set frame to use a world axis as a tangent and
- return*/
- if (!self->parent || dot >= 1.0-DBL_EPSILON*3) {
+ if (!self->parent || dot >= 1.0-DBL_EPSILON*2) {
self->head.result = 0;
return;
}
- mul_v3_d(vec2, l2);
- add_v3_v3_d(vec2, v1->tx);
- solve_bending_gradient(v2->txconst, v1->tx, vec2, g1, g2, g3, 0.0);
+ mul_v3_d(vec1, l2);
+ add_v3_v3_d(vec1, v1->tx);
+ solve_bending_gradient(v2->tx, v1->tx, vec1, g1, g2, g3, 0.0);
+
if (isnan(g1[0])) {
self->head.result = 0.0;
return;
}
- copy_v3_v3_d(self->head.grad[1], g2);
- copy_v3_v3_d(self->head.grad[0], g3);
+ mul_v3_d(g1, -1.0);
+ copy_v3_v3_d(self->head.grad[0], g1);
- if (len_v3_d(self->head.grad[0]) > 3.0 || len_v3_d(self->head.grad[1]) > 3.0) {
- //printf("eek!\n");
- }
- self->head.result = fabs(saacos_d(dot));
+ //if (len_v3_d(self->head.grad[0]) > 3.0 || len_v3_d(self->head.grad[1]) > 3.0) {
+ //printf("eek!\n");
+ //}
+ self->head.result = -saacos(dot); //dot < 0.0 ? -saacos_d(dot) : saacos_d(dot);
return;
}
@@ -389,14 +416,14 @@
static GoalConstraint *make_goal_constraint(ClothVertex *v1, ClothVertex *v2,
- GoalConstraint *last, float k, float restlen, MemArena *arena)
+ GoalConstraint *last, float k, float restlen, MemArena *arena)
{
GoalConstraint *con = BLI_memarena_alloc(arena, sizeof(*con));
memset(con, 0, sizeof(*con));
con->head.verts[0] = v2;
- con->head.verts[1] = v1;
+ //con->head.verts[1] = v1;
//if (last)
// con->head.verts[2] = last->head.verts[1];
@@ -416,38 +443,35 @@
typedef struct BendingConstraint {
Constraint head;
+ ClothVertex *v2, *v3;
+
cdouble_t restlen;
} BendingConstraint;
static void eval_bending_constraint(void *vself)
{
BendingConstraint *self = vself;
- ClothVertex *v1 = self->head.verts[0], *v2=self->head.verts[1], *v3=self->head.verts[2];
- cdouble_t quat[4], vec1[3], vec2[3], tan[3], vec5[3], nor[3], dot, l2;
- cdouble_t g1[3], g2[3], g3[3], t, rest;
+ ClothVertex *v1 = self->head.verts[0], *v2=self->v2, *v3=self->v3;
+ cdouble_t vec1[3], vec2[3], dot, l2;
+ cdouble_t g1[3], g2[3], g3[3], t;
- v3 = self->head.verts[2];
-
/*calculate reference frame*/
sub_v3_v3v3_d(vec1, v1->tx, v2->tx);
sub_v3_v3v3_d(vec2, v3->tx, v2->tx);
- l2 = len_v3_d(vec1);
- normalize_v3_d(vec1);
- normalize_v3_d(vec2);
+ l2 = dot_v3_d(vec1);
+ normalize_fast_d(vec1);
+ normalize_fast_d(vec2);
+
dot = dot_v3v3_d(vec1, vec2);
t = saacos_d(dot);
- if (t > self->restlen || fabs(t - self->restlen) < 0.001 || dot <= -0.99999) {
+ if (t > self->restlen || fabs(t - self->restlen) < 0.0001 || dot <= -0.99999) {
self->head.result = 0.0;
return;
}
- memset(g1, 0, sizeof(g1));
- memset(g2, 0, sizeof(g2));
- memset(g3, 0, sizeof(g3));
-
solve_bending_gradient(v1->tx, v2->tx, v3->tx, g1, g2, g3, self->restlen);
if (isnan(g1[0])) {
@@ -455,14 +479,18 @@
return;
}
+ //mul_v3_d(g1, 1.0/3.0);
+ //mul_v3_d(g2, 1.0/3.0);
+ //mul_v3_d(g3, 1.0/3.0);
+
copy_v3_v3_d(self->head.grad[0], g1);
copy_v3_v3_d(self->head.grad[1], g2);
copy_v3_v3_d(self->head.grad[2], g3);
- self->head.result = -fabs(t - self->restlen);
+ self->head.result = t - self->restlen;
- if (fabs(self->head.result) < 0.00001)
- self->head.result = 0.0;
+ //if (fabs(self->head.result) < 0.00001)
+ // self->head.result = 0.0;
return;
}
@@ -477,16 +505,19 @@
con->head.verts[0] = v1;
con->head.verts[1] = v2;
- con->head.verts[2] = v3;
+ //con->head.verts[2] = v3;
+ con->v2 = v2;
+ con->v3 = v3;
+
con->head.eval = eval_bending_constraint;
con->head.type = CON_BEND;
- sub_v3_v3v3_d(vec1, v1->txconst, v2->txconst);
- sub_v3_v3v3_d(vec2, v3->txconst, v2->txconst);
+ sub_v3_v3v3_d(vec1, v1->xrest, v2->xrest);
+ sub_v3_v3v3_d(vec2, v3->xrest, v2->xrest);
- normalize_v3_d(vec1);
- normalize_v3_d(vec2);
+ normalize_fast_d(vec1);
+ normalize_fast_d(vec2);
con->restlen = saacos_d(dot_v3v3_d(vec1, vec2));
con->head.k = k;
@@ -749,11 +780,11 @@
VECSUB(off, cent, ocent);
if (off[0] != 0.0 || off[1] != 0.0 || off[2] != 0.0) {
- normalize_v3_d(off);
+ normalize_fast_d(off);
VECMULS(off, MAX2(thickness, 0.01));
} else {
VECCOPY(off, no)
- VECMULS(off, MAX2(thickness, 0.01));
+ VECMULS(off, MAX2(thickness, 0.01));
}
eps = 1.0001; // + DBL_EPSILON*2;
@@ -793,8 +824,8 @@
/*expands faces by eps*/
static int swept_tri_overlap(cdouble_t tverts1[3][3], cdouble_t tlastverts1[3][3],
- cdouble_t tverts2[3][3], cdouble_t tlastverts2[3][3],
- cdouble_t eps, cdouble_t thickness)
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list