More F-Modifier Tweaks:
[blender.git] / source / blender / blenkernel / intern / implicit.c
index fc5213d5532faca3508f0f3cf08281317d0c8c44..6912a65886a076cacbafd65ae5b14c988f882d51 100644 (file)
@@ -15,7 +15,7 @@
 *
 * 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * The Original Code is Copyright (C) Blender Foundation
 * All rights reserved.
@@ -33,6 +33,7 @@
 
 #include "DNA_cloth_types.h"   
 #include "DNA_scene_types.h"
+#include "DNA_object_force.h"
 
 #include "BKE_effect.h"
 #include "BKE_global.h"
@@ -43,7 +44,7 @@
 #include <windows.h>
 static LARGE_INTEGER _itstart, _itend;
 static LARGE_INTEGER ifreq;
-void itstart(void)
+static void itstart(void)
 {
        static int first = 1;
        if(first) {
@@ -52,7 +53,7 @@ void itstart(void)
        }
        QueryPerformanceCounter(&_itstart);
 }
-void itend(void)
+static void itend(void)
 {
        QueryPerformanceCounter(&_itend);
 }
@@ -74,7 +75,7 @@ double itval()
 {
        gettimeofday(&_itstart, &itz);
 }
-void itend(void)
+static void itend(void)
 {
        gettimeofday(&_itend,&itz);
 }
@@ -155,7 +156,7 @@ DO_INLINE void mul_fvectorT_fvectorS(float to[3][3], float vectorA[3], float vec
 
 
 /* printf vector[3] on console: for debug output */
-void print_fvector(float m3[3])
+static void print_fvector(float m3[3])
 {
        printf("%f\n%f\n%f\n\n",m3[0],m3[1],m3[2]);
 }
@@ -296,13 +297,15 @@ DO_INLINE void sub_lfvector_lfvector(float (*to)[3], float (*fLongVectorA)[3], f
 ///////////////////////////
 // 3x3 matrix
 ///////////////////////////
+#if 0
 /* printf 3x3 matrix on console: for debug output */
-void print_fmatrix(float m3[3][3])
+static void print_fmatrix(float m3[3][3])
 {
        printf("%f\t%f\t%f\n",m3[0][0],m3[0][1],m3[0][2]);
        printf("%f\t%f\t%f\n",m3[1][0],m3[1][1],m3[1][2]);
        printf("%f\t%f\t%f\n\n",m3[2][0],m3[2][1],m3[2][2]);
 }
+#endif
 
 /* copy 3x3 matrix */
 DO_INLINE void cp_fmatrix(float to[3][3], float from[3][3])
@@ -496,7 +499,8 @@ DO_INLINE void mulsub_fmatrix_fvector(float to[3], float matrix[3][3], float fro
 // SPARSE SYMMETRIC big matrix with 3x3 matrix entries
 ///////////////////////////
 /* printf a big matrix on console: for debug output */
-void print_bfmatrix(fmatrix3x3 *m3)
+#if 0
+static void print_bfmatrix(fmatrix3x3 *m3)
 {
        unsigned int i = 0;
 
@@ -505,6 +509,8 @@ void print_bfmatrix(fmatrix3x3 *m3)
                print_fmatrix(m3[i].m);
        }
 }
+#endif
+
 /* create big matrix */
 DO_INLINE fmatrix3x3 *create_bfmatrix(unsigned int verts, unsigned int springs)
 {
@@ -887,7 +893,7 @@ DO_INLINE void filter(lfVector *V, fmatrix3x3 *S)
        }
 }
 
-int  cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S)
+static int  cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S)
 {
        // Solves for unknown X in equation AX=B
        unsigned int conjgrad_loopcount=0, conjgrad_looplimit=100;
@@ -968,9 +974,10 @@ DO_INLINE void BuildPPinv(fmatrix3x3 *lA, fmatrix3x3 *P, fmatrix3x3 *Pinv)
                
        }
 }
+#if 0
 /*
 // version 1.3
-int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, fmatrix3x3 *P, fmatrix3x3 *Pinv)
+static int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, fmatrix3x3 *P, fmatrix3x3 *Pinv)
 {
        unsigned int numverts = lA[0].vcount, iterations = 0, conjgrad_looplimit=100;
        float delta0 = 0, deltaNew = 0, deltaOld = 0, alpha = 0;
@@ -1038,7 +1045,7 @@ int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fma
 }
 */
 // version 1.4
-int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, fmatrix3x3 *P, fmatrix3x3 *Pinv, fmatrix3x3 *bigI)
+static int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, fmatrix3x3 *P, fmatrix3x3 *Pinv, fmatrix3x3 *bigI)
 {
        unsigned int numverts = lA[0].vcount, iterations = 0, conjgrad_looplimit=100;
        float delta0 = 0, deltaNew = 0, deltaOld = 0, alpha = 0, tol = 0;
@@ -1139,6 +1146,7 @@ int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fma
        
        return iterations<conjgrad_looplimit;
 }
+#endif
 
 // outer product is NOT cross product!!!
 DO_INLINE void dfdx_spring_type1(float to[3][3], float extent[3], float length, float L, float dot, float k)
@@ -1183,7 +1191,8 @@ DO_INLINE void dfdx_spring(float to[3][3],  float dir[3],float length,float L,fl
        //return  ( (I-outerprod(dir,dir))*Min(1.0f,rest/length) - I) * -k;
        mul_fvectorT_fvector(to, dir, dir);
        sub_fmatrix_fmatrix(to, I, to);
-       mul_fmatrix_S(to, (((L/length)> 1.0f) ? (1.0f): (L/length))); 
+
+       mul_fmatrix_S(to, (L/length)); 
        sub_fmatrix_fmatrix(to, to, I);
        mul_fmatrix_S(to, -k);
 }
@@ -1218,6 +1227,8 @@ DO_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
        float nulldfdx[3][3]={ {0,0,0}, {0,0,0}, {0,0,0}};
        
        float scaling = 0.0;
+
+       int no_compress = clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_NO_SPRING_COMPRESS;
        
        VECCOPY(s->f, nullf);
        cp_fmatrix(s->dfdx, nulldfdx);
@@ -1254,7 +1265,7 @@ DO_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
        // calculate force of structural + shear springs
        if((s->type & CLOTH_SPRING_TYPE_STRUCTURAL) || (s->type & CLOTH_SPRING_TYPE_SHEAR))
        {
-               if(length > L) // only on elonglation
+               if(length > L || no_compress)
                {
                        s->flags |= CLOTH_SPRING_FLAG_NEEDED;
                        
@@ -1388,26 +1399,166 @@ static void CalcFloat4( float *v1, float *v2, float *v3, float *v4, float *n)
        n[2]= n1[0]*n2[1]-n1[1]*n2[0];
 }
 
-float calculateVertexWindForce(float wind[3], float vertexnormal[3])  
+static float calculateVertexWindForce(float wind[3], float vertexnormal[3])  
 {
        return (INPR(wind, vertexnormal));
 }
 
+typedef struct HairGridVert {
+       float velocity[3];
+       float density;
+} HairGridVert;
+#define HAIR_GRID_INDEX(vec, min, max, axis) (int)( (vec[axis] - min[axis]) / (max[axis] - min[axis]) * 9.99f );
+/* Smoothing of hair velocities:
+ * adapted from
+               Volumetric Methods for Simulation and Rendering of Hair
+               by Lena Petrovic, Mark Henne and John Anderson
+ *             Pixar Technical Memo #06-08, Pixar Animation Studios
+ */
+static void hair_velocity_smoothing(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVector *lV, int numverts)
+{
+       /* TODO: This is an initial implementation and should be made much better in due time.
+        * What should at least be implemented is a grid size parameter and a smoothing kernel
+        * for bigger grids.
+        */
+
+       /* 10x10x10 grid gives nice initial results */
+       HairGridVert grid[10][10][10];
+       HairGridVert colg[10][10][10];
+       ListBase *colliders = get_collider_cache(clmd->scene, NULL);
+       ColliderCache *col = NULL;
+       float gmin[3], gmax[3], density;
+       /* 2.0f is an experimental value that seems to give good results */
+       float smoothfac = 2.0f * clmd->sim_parms->velocity_smooth;
+       float collfac = 2.0f * clmd->sim_parms->collider_friction;
+       int     v = 0;
+       int     i = 0;
+       int     j = 0;
+       int     k = 0;
+
+       INIT_MINMAX(gmin, gmax);
+
+       for(i = 0; i < numverts; i++)
+               DO_MINMAX(lX[i], gmin, gmax);
+
+       /* initialize grid */
+       for(i = 0; i < 10; i++) {
+               for(j = 0; j < 10; j++) {
+                       for(k = 0; k < 10; k++) {
+                               grid[i][j][k].velocity[0] = 0.0f;
+                               grid[i][j][k].velocity[1] = 0.0f;
+                               grid[i][j][k].velocity[2] = 0.0f;
+                               grid[i][j][k].density = 0.0f;
+
+                               colg[i][j][k].velocity[0] = 0.0f;
+                               colg[i][j][k].velocity[1] = 0.0f;
+                               colg[i][j][k].velocity[2] = 0.0f;
+                               colg[i][j][k].density = 0.0f;
+                       }
+               }
+       }
+
+       /* gather velocities & density */
+       if(smoothfac > 0.0f) for(v = 0; v < numverts; v++) {
+               i = HAIR_GRID_INDEX(lX[v], gmin, gmax, 0);
+               j = HAIR_GRID_INDEX(lX[v], gmin, gmax, 1);
+               k = HAIR_GRID_INDEX(lX[v], gmin, gmax, 2);
+
+               grid[i][j][k].velocity[0] += lV[v][0];
+               grid[i][j][k].velocity[1] += lV[v][1];
+               grid[i][j][k].velocity[2] += lV[v][2];
+               grid[i][j][k].density += 1.0f;
+       }
+
+       /* gather colliders */
+       if(colliders && collfac > 0.0f) for(col = colliders->first; col; col = col->next)
+       {
+               MVert *loc0 = col->collmd->x;
+               MVert *loc1 = col->collmd->xnew;
+               float vel[3];
+
+               for(v=0; v<col->collmd->numverts; v++, loc0++, loc1++) {
+                       i = HAIR_GRID_INDEX(loc1->co, gmin, gmax, 0);
+
+                       if(i>=0 && i<10) {
+                               j = HAIR_GRID_INDEX(loc1->co, gmin, gmax, 1);
+
+                               if(j>=0 && j<10) {
+                                       k = HAIR_GRID_INDEX(loc1->co, gmin, gmax, 2);
+
+                                       if(k>=0 && k<10) {
+                                               VECSUB(vel, loc1->co, loc0->co);
+
+                                               colg[i][j][k].velocity[0] += vel[0];
+                                               colg[i][j][k].velocity[1] += vel[1];
+                                               colg[i][j][k].velocity[2] += vel[2];
+                                               colg[i][j][k].density += 1.0;
+                                       }
+                               }
+                       }
+               }
+       }
+       
+
+       /* divide velocity with density */
+       for(i = 0; i < 10; i++) {
+               for(j = 0; j < 10; j++) {
+                       for(k = 0; k < 10; k++) {
+                               density = grid[i][j][k].density;
+                               if(density > 0.0f) {
+                                       grid[i][j][k].velocity[0] /= density;
+                                       grid[i][j][k].velocity[1] /= density;
+                                       grid[i][j][k].velocity[2] /= density;
+                               }
+
+                               density = colg[i][j][k].density;
+                               if(density > 0.0f) {
+                                       colg[i][j][k].velocity[0] /= density;
+                                       colg[i][j][k].velocity[1] /= density;
+                                       colg[i][j][k].velocity[2] /= density;
+                               }
+                       }
+               }
+       }
+
+       /* calculate forces */
+       for(v = 0; v < numverts; v++) {
+               i = HAIR_GRID_INDEX(lX[v], gmin, gmax, 0);
+               j = HAIR_GRID_INDEX(lX[v], gmin, gmax, 1);
+               k = HAIR_GRID_INDEX(lX[v], gmin, gmax, 2);
+
+               lF[v][0] += smoothfac * (grid[i][j][k].velocity[0] - lV[v][0]);
+               lF[v][1] += smoothfac * (grid[i][j][k].velocity[1] - lV[v][1]);
+               lF[v][2] += smoothfac * (grid[i][j][k].velocity[2] - lV[v][2]);
+
+               if(colg[i][j][k].density > 0.0f) {
+                       lF[v][0] += collfac * (colg[i][j][k].velocity[0] - lV[v][0]);
+                       lF[v][1] += collfac * (colg[i][j][k].velocity[1] - lV[v][1]);
+                       lF[v][2] += collfac * (colg[i][j][k].velocity[2] - lV[v][2]);
+               }
+       }
+
+       free_collider_cache(&colliders);
+}
 static void cloth_calc_force(ClothModifierData *clmd, float frame, lfVector *lF, lfVector *lX, lfVector *lV, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, ListBase *effectors, float time, fmatrix3x3 *M)
 {
        /* Collect forces and derivatives:  F,dFdX,dFdV */
        Cloth           *cloth          = clmd->clothObject;
        int             i               = 0;
        float           spring_air      = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
-       float           gravity[3];
+       float           gravity[3] = {0.0f, 0.0f, 0.0f};
        float           tm2[3][3]       = {{-spring_air,0,0}, {0,-spring_air,0},{0,0,-spring_air}};
        MFace           *mfaces         = cloth->mfaces;
        unsigned int numverts = cloth->numverts;
        LinkNode *search = cloth->springs;
        lfVector *winvec;
+       EffectedPoint epoint;
 
-       VECCOPY(gravity, clmd->sim_parms->gravity);
-       mul_fvector_S(gravity, gravity, 0.001f); /* scale gravity force */
+       /* global acceleration (gravitation) */
+       if(clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) {
+               VECCOPY(gravity, clmd->scene->physics_settings.gravity);
+               mul_fvector_S(gravity, gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity); /* scale gravity force */
+       }
 
        /* set dFdX jacobi matrix to zero */
        init_bfmatrix(dFdX, ZERO);
@@ -1416,6 +1567,9 @@ static void cloth_calc_force(ClothModifierData *clmd, float frame, lfVector *lF,
 
        init_lfvector(lF, gravity, numverts);
        
+       if(clmd->sim_parms->velocity_smooth > 0.0f || clmd->sim_parms->collider_friction > 0.0f)
+               hair_velocity_smoothing(clmd, lF, lX, lV, numverts);
+
        /* multiply lF with mass matrix
        // force = mass * acceleration (in this case: gravity)
        */
@@ -1439,10 +1593,9 @@ static void cloth_calc_force(ClothModifierData *clmd, float frame, lfVector *lF,
                
                // precalculate wind forces
                for(i = 0; i < cloth->numverts; i++)
-               {
-                       float speed[3] = {0.0f, 0.0f,0.0f};
-                       
-                       pdDoEffectors(clmd->scene, effectors, lX[i], winvec[i], speed, frame, 0.0f, 0);
+               {       
+                       pd_point_from_loc(clmd->scene, (float*)lX[i], (float*)lV[i], i, &epoint);
+                       pdDoEffectors(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL);
                }
                
                for(i = 0; i < cloth->numfaces; i++)
@@ -1460,31 +1613,61 @@ static void cloth_calc_force(ClothModifierData *clmd, float frame, lfVector *lF,
                                CalcFloat(lX[mfaces[i].v1],lX[mfaces[i].v2],lX[mfaces[i].v3],triunnormal);
                        
                        VECCOPY(trinormal, triunnormal);
-                       Normalize(trinormal);
+                       normalize_v3(trinormal);
                        
                        // add wind from v1
                        VECCOPY(tmp, trinormal);
-                       VecMulf(tmp, calculateVertexWindForce(winvec[mfaces[i].v1], triunnormal));
+                       mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v1], triunnormal));
                        VECADDS(lF[mfaces[i].v1], lF[mfaces[i].v1], tmp, factor);
                        
                        // add wind from v2
                        VECCOPY(tmp, trinormal);
-                       VecMulf(tmp, calculateVertexWindForce(winvec[mfaces[i].v2], triunnormal));
+                       mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v2], triunnormal));
                        VECADDS(lF[mfaces[i].v2], lF[mfaces[i].v2], tmp, factor);
                        
                        // add wind from v3
                        VECCOPY(tmp, trinormal);
-                       VecMulf(tmp, calculateVertexWindForce(winvec[mfaces[i].v3], triunnormal));
+                       mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v3], triunnormal));
                        VECADDS(lF[mfaces[i].v3], lF[mfaces[i].v3], tmp, factor);
                        
                        // add wind from v4
                        if(mfaces[i].v4)
                        {
                                VECCOPY(tmp, trinormal);
-                               VecMulf(tmp, calculateVertexWindForce(winvec[mfaces[i].v4], triunnormal));
+                               mul_v3_fl(tmp, calculateVertexWindForce(winvec[mfaces[i].v4], triunnormal));
                                VECADDS(lF[mfaces[i].v4], lF[mfaces[i].v4], tmp, factor);
                        }
                }
+
+               /* Hair has only edges */
+               if(cloth->numfaces == 0) {
+                       ClothSpring *spring;
+                       float edgevec[3]={0,0,0}; //edge vector
+                       float edgeunnormal[3]={0,0,0}; // not-normalized-edge normal
+                       float tmp[3]={0,0,0};
+                       float factor = 0.01;
+
+                       search = cloth->springs;
+                       while(search) {
+                               spring = search->link;
+                               
+                               if(spring->type == CLOTH_SPRING_TYPE_STRUCTURAL) {
+                                       VECSUB(edgevec, (float*)lX[spring->ij], (float*)lX[spring->kl]);
+
+                                       project_v3_v3v3(tmp, winvec[spring->ij], edgevec);
+                                       VECSUB(edgeunnormal, winvec[spring->ij], tmp);
+                                       /* hair doesn't stretch too much so we can use restlen pretty safely */
+                                       VECADDS(lF[spring->ij], lF[spring->ij], edgeunnormal, spring->restlen * factor);
+
+                                       project_v3_v3v3(tmp, winvec[spring->kl], edgevec);
+                                       VECSUB(edgeunnormal, winvec[spring->kl], tmp);
+                                       VECADDS(lF[spring->kl], lF[spring->kl], edgeunnormal, spring->restlen * factor);
+                               }
+
+                               search = search->next;
+                       }
+               }
+
                del_lfvector(winvec);
        }
                
@@ -1511,7 +1694,7 @@ static void cloth_calc_force(ClothModifierData *clmd, float frame, lfVector *lF,
        // printf("\n");
 }
 
-void simulate_implicit_euler(lfVector *Vnew, lfVector *lX, lfVector *lV, lfVector *lF, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, float dt, fmatrix3x3 *A, lfVector *B, lfVector *dV, fmatrix3x3 *S, lfVector *z, lfVector *olddV, fmatrix3x3 *P, fmatrix3x3 *Pinv, fmatrix3x3 *M, fmatrix3x3 *bigI)
+static void simulate_implicit_euler(lfVector *Vnew, lfVector *lX, lfVector *lV, lfVector *lF, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, float dt, fmatrix3x3 *A, lfVector *B, lfVector *dV, fmatrix3x3 *S, lfVector *z, lfVector *olddV, fmatrix3x3 *P, fmatrix3x3 *Pinv, fmatrix3x3 *M, fmatrix3x3 *bigI)
 {
        unsigned int numverts = dFdV[0].vcount;
 
@@ -1562,7 +1745,7 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                        if(verts [i].flags & CLOTH_VERT_FLAG_PINNED)
                        {                       
                                VECSUB(id->V[i], verts[i].xconst, verts[i].xold);
-                               // VecMulf(id->V[i], clmd->sim_parms->stepsPerFrame);
+                               // mul_v3_fl(id->V[i], clmd->sim_parms->stepsPerFrame);
                        }
                }       
        }
@@ -1570,9 +1753,7 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
        while(step < tf)
        {       
                // calculate forces
-               effectors= pdInitEffectors(clmd->scene, ob, NULL);
                cloth_calc_force(clmd, frame, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step, id->M);
-               if(effectors) pdEndEffectors(effectors);
                
                // calculate new velocity
                simulate_implicit_euler(id->Vnew, id->X, id->V, id->F, id->dFdV, id->dFdX, dt, id->A, id->B, id->dV, id->S, id->z, id->olddV, id->P, id->Pinv, id->M, id->bigI);
@@ -1597,8 +1778,8 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                        
                        VECCOPY(verts[i].txold, id->X[i]);
                }
-               
-               if(clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED)
+
+               if(clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED && clmd->clothObject->bvhtree)
                {
                        float temp = clmd->sim_parms->stepsPerFrame;
                        /* not too nice hack, but collisions need this correction -jahka */
@@ -1637,7 +1818,7 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                                        
                                        VECCOPY(id->Xnew[i], verts[i].tx);
                                        VECCOPY(id->Vnew[i], verts[i].tv);
-                                       VecMulf(id->Vnew[i], clmd->sim_parms->stepsPerFrame);
+                                       mul_v3_fl(id->Vnew[i], clmd->sim_parms->stepsPerFrame);
                                }
                        }
                        
@@ -1655,9 +1836,7 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                                cp_lfvector(id->V, id->Vnew, numverts);
                                
                                // calculate 
-                               effectors= pdInitEffectors(clmd->scene, ob, NULL);
                                cloth_calc_force(clmd, frame, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step+dt, id->M);      
-                               if(effectors) pdEndEffectors(effectors);
                                
                                simulate_implicit_euler(id->Vnew, id->X, id->V, id->F, id->dFdV, id->dFdX, dt / 2.0f, id->A, id->B, id->dV, id->S, id->z, id->olddV, id->P, id->Pinv, id->M, id->bigI);
                        }