More F-Modifier Tweaks:
[blender.git] / source / blender / blenkernel / intern / implicit.c
index 96c3c33afb86e3d16a705b3777d990bea05bf0d6..6912a65886a076cacbafd65ae5b14c988f882d51 100644 (file)
@@ -1,15 +1,12 @@
 /*  implicit.c      
 * 
 *
-* ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
+* ***** 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. The Blender
-* Foundation also sells licenses for use in proprietary software under
-* the Blender License.  See http://www.blender.org/BL/ for information
-* about this.
+* 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
@@ -18,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.
 *
 * Contributor(s): none yet.
 *
-* ***** END GPL/BL DUAL LICENSE BLOCK *****
+* ***** END GPL LICENSE BLOCK *****
 */
-#include "math.h"
-#include "float.h"
-#include <stdlib.h>
-#include <string.h>
-#include <stdio.h>
+
 #include "MEM_guardedalloc.h"
-/* types */
-#include "DNA_curve_types.h"
-#include "DNA_object_types.h"
-#include "DNA_object_force.h"  
+
+#include "BKE_cloth.h"
+
 #include "DNA_cloth_types.h"   
-#include "DNA_key_types.h"
-#include "DNA_mesh_types.h"
-#include "DNA_modifier_types.h"
-#include "DNA_meshdata_types.h"
-#include "DNA_lattice_types.h"
 #include "DNA_scene_types.h"
-#include "DNA_modifier_types.h"
-#include "BLI_blenlib.h"
-#include "BLI_arithb.h"
-#include "BLI_threads.h"
-#include "BKE_curve.h"
-#include "BKE_displist.h"
+#include "DNA_object_force.h"
+
 #include "BKE_effect.h"
 #include "BKE_global.h"
-#include "BKE_key.h"
-#include "BKE_object.h"
 #include "BKE_cloth.h"
-#include "BKE_modifier.h"
 #include "BKE_utildefines.h"
-#include "BKE_global.h"
-#include "BIF_editdeform.h"
-
 
 #ifdef _WIN32
 #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) {
@@ -76,7 +53,7 @@ void itstart(void)
        }
        QueryPerformanceCounter(&_itstart);
 }
-void itend(void)
+static void itend(void)
 {
        QueryPerformanceCounter(&_itend);
 }
@@ -98,7 +75,7 @@ double itval()
 {
        gettimeofday(&_itstart, &itz);
 }
-void itend(void)
+static void itend(void)
 {
        gettimeofday(&_itend,&itz);
 }
@@ -179,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]);
 }
@@ -320,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])
@@ -520,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;
 
@@ -529,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)
 {
@@ -911,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;
@@ -992,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;
@@ -1062,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;
@@ -1163,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)
@@ -1207,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);
 }
@@ -1242,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);
@@ -1278,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;
                        
@@ -1378,63 +1365,200 @@ DO_INLINE void cloth_apply_spring_force(ClothModifierData *clmd, ClothSpring *s,
        }       
 }
 
-DO_INLINE void calculateTriangleNormal(float to[3], lfVector *X, MFace mface)
+
+static void CalcFloat( float *v1, float *v2, float *v3, float *n)
 {
-       float v1[3], v2[3];
+       float n1[3],n2[3];
 
-       VECSUB(v1, X[mface.v2], X[mface.v1]);
-       VECSUB(v2, X[mface.v3], X[mface.v1]);
-       cross_fvector(to, v1, v2);
+       n1[0]= v1[0]-v2[0];
+       n2[0]= v2[0]-v3[0];
+       n1[1]= v1[1]-v2[1];
+       n2[1]= v2[1]-v3[1];
+       n1[2]= v1[2]-v2[2];
+       n2[2]= v2[2]-v3[2];
+       n[0]= n1[1]*n2[2]-n1[2]*n2[1];
+       n[1]= n1[2]*n2[0]-n1[0]*n2[2];
+       n[2]= n1[0]*n2[1]-n1[1]*n2[0];
 }
 
-DO_INLINE void calculatQuadNormal(float to[3], lfVector *X, MFace mface)
+static void CalcFloat4( float *v1, float *v2, float *v3, float *v4, float *n)
 {
-       float temp = CalcNormFloat4(X[mface.v1],X[mface.v2],X[mface.v3],X[mface.v4],to);
-       mul_fvector_S(to, to, temp);
+       /* real cross! */
+       float n1[3],n2[3];
+
+       n1[0]= v1[0]-v3[0];
+       n1[1]= v1[1]-v3[1];
+       n1[2]= v1[2]-v3[2];
+
+       n2[0]= v2[0]-v4[0];
+       n2[1]= v2[1]-v4[1];
+       n2[2]= v2[2]-v4[2];
+
+       n[0]= n1[1]*n2[2]-n1[2]*n2[1];
+       n[1]= n1[2]*n2[0]-n1[0]*n2[2];
+       n[2]= n1[0]*n2[1]-n1[1]*n2[0];
 }
 
-void calculateWeightedVertexNormal(ClothModifierData *clmd, MFace *mfaces, float to[3], int index, lfVector *X)
+static float calculateVertexWindForce(float wind[3], float vertexnormal[3])  
 {
-       float temp[3]; 
-       int i;
-       Cloth *cloth = clmd->clothObject;
+       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);
 
-       for(i = 0; i < cloth->numfaces; i++)
+               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)
        {
-               // check if this triangle contains the selected vertex
-               if(mfaces[i].v1 == index || mfaces[i].v2 == index || mfaces[i].v3 == index || mfaces[i].v4 == index)
-               {
-                       calculatQuadNormal(temp, X, mfaces[i]);
-                       VECADD(to, to, temp);
+               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;
+                                       }
+                               }
+                       }
                }
        }
-}
-float calculateVertexWindForce(float wind[3], float vertexnormal[3])  
-{
-       return fabs(INPR(wind, vertexnormal) * 0.5f);
-}
+       
 
-DO_INLINE void calc_triangle_force(ClothModifierData *clmd, MFace mface, lfVector *F, lfVector *X, lfVector *V, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, ListBase *effectors)
-{      
+       /* 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;
+                               }
+                       }
+               }
+       }
 
-void cloth_calc_force(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVector *lV, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, ListBase *effectors, float time, fmatrix3x3 *M)
+       /* 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;
-       long            i               = 0;
+       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;
-       float wind_normalized[3];
        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);
@@ -1443,8 +1567,13 @@ void cloth_calc_force(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVec
 
        init_lfvector(lF, gravity, numverts);
        
-       // multiply lF with mass matrix
-       for(i = 0; i < (long)numverts; i++)
+       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)
+       */
+       for(i = 0; i < numverts; i++)
        {
                float temp[3];
                VECCOPY(temp, lF[i]);
@@ -1455,27 +1584,91 @@ void cloth_calc_force(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVec
        
        /* handle external forces like wind */
        if(effectors)
-       {
-               float speed[3] = {0.0f, 0.0f,0.0f};
-               float force[3]= {0.0f, 0.0f, 0.0f};
+       {       
+               // 0 = force, 1 = normalized force
+               winvec = create_lfvector(cloth->numverts);
                
-#pragma omp parallel for private (i) shared(lF)
-               for(i = 0; i < (long)(cloth->numverts); i++)
+               if(!winvec)
+                       printf("winvec: out of memory in implicit.c\n");
+               
+               // precalculate wind forces
+               for(i = 0; i < cloth->numverts; i++)
+               {       
+                       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++)
                {
-                       float vertexnormal[3]={0,0,0};
-                       float fieldfactor = 1000.0f; // windfactor  = 250.0f; // from sb
+                       float trinormal[3]={0,0,0}; // normalized triangle normal
+                       float triunnormal[3]={0,0,0}; // not-normalized-triangle normal
+                       float tmp[3]={0,0,0};
+                       float factor = (mfaces[i].v4) ? 0.25 : 1.0 / 3.0;
+                       factor *= 0.02;
                        
-                       pdDoEffectors(effectors, lX[i], force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);          
+                       // calculate face normal
+                       if(mfaces[i].v4)
+                               CalcFloat4(lX[mfaces[i].v1],lX[mfaces[i].v2],lX[mfaces[i].v3],lX[mfaces[i].v4],triunnormal);
+                       else
+                               CalcFloat(lX[mfaces[i].v1],lX[mfaces[i].v2],lX[mfaces[i].v3],triunnormal);
                        
-                       // TODO apply forcefields here
-                       VECADDS(lF[i], lF[i], force, fieldfactor*0.01f);
-
-                       VECCOPY(wind_normalized, speed);
-                       Normalize(wind_normalized);
+                       VECCOPY(trinormal, triunnormal);
+                       normalize_v3(trinormal);
                        
-                       calculateWeightedVertexNormal(clmd, mfaces, vertexnormal, i, lX);
-                       VECADDS(lF[i], lF[i], wind_normalized, -calculateVertexWindForce(speed, vertexnormal));
+                       // add wind from v1
+                       VECCOPY(tmp, trinormal);
+                       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);
+                       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);
+                       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);
+                               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);
        }
                
        // calculate spring forces
@@ -1501,7 +1694,7 @@ void cloth_calc_force(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVec
        // 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;
 
@@ -1536,11 +1729,11 @@ void simulate_implicit_euler(lfVector *Vnew, lfVector *lX, lfVector *lV, lfVecto
 int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors)
 {              
        unsigned int i=0;
-       float step=0.0f, tf=1.0f;
+       float step=0.0f, tf=clmd->sim_parms->timescale;
        Cloth *cloth = clmd->clothObject;
        ClothVertex *verts = cloth->verts;
        unsigned int numverts = cloth->numverts;
-       float dt = 1.0f / clmd->sim_parms->stepsPerFrame;
+       float dt = clmd->sim_parms->timescale / clmd->sim_parms->stepsPerFrame;
        Implicit_Data *id = cloth->implicit;
        int result = 0;
        
@@ -1552,17 +1745,15 @@ 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);
                        }
                }       
        }
        
        while(step < tf)
        {       
-               effectors= pdInitEffectors(ob,NULL);
-               
                // calculate forces
-               cloth_calc_force(clmd, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step, id->M);
+               cloth_calc_force(clmd, frame, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step, id->M);
                
                // 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);
@@ -1587,9 +1778,13 @@ 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 */
+                       clmd->sim_parms->stepsPerFrame /= clmd->sim_parms->timescale;
+
                        // collisions 
                        // itstart();
                        
@@ -1601,10 +1796,17 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                                VECSUB(verts[i].tv, verts[i].tx, verts[i].txold);
                                VECCOPY(verts[i].v, verts[i].tv);
                        }
-       
+                       
                        // call collision function
-                       result = cloth_bvh_objcollision(clmd, step + dt, dt);
-       
+                       // TODO: check if "step" or "step+dt" is correct - dg
+                       result = cloth_bvh_objcollision(ob, clmd, step/clmd->sim_parms->timescale, dt/clmd->sim_parms->timescale);
+                       
+                       // correct velocity again, just to be sure we had to change it due to adaptive collisions
+                       for(i = 0; i < numverts; i++)
+                       {
+                               VECSUB(verts[i].tv, verts[i].tx, id->X[i]);
+                       }
+                       
                        // copy corrected positions back to simulation
                        for(i = 0; i < numverts; i++)
                        {               
@@ -1616,10 +1818,13 @@ 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);
                                }
                        }
                        
+                       /* restore original stepsPerFrame */
+                       clmd->sim_parms->stepsPerFrame = temp;
+                       
                        // X = Xnew;
                        cp_lfvector(id->X, id->Xnew, numverts);
                        
@@ -1631,10 +1836,10 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                                cp_lfvector(id->V, id->Vnew, numverts);
                                
                                // calculate 
-                               cloth_calc_force(clmd, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step+dt, id->M);     
+                               cloth_calc_force(clmd, frame, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step+dt, id->M);      
+                               
                                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);
                        }
-                       
                }
                else
                {
@@ -1649,8 +1854,6 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                cp_lfvector(id->V, id->Vnew, numverts);
                
                step += dt;
-
-               if(effectors) pdEndEffectors(effectors);
                
        }
 
@@ -1669,6 +1872,7 @@ int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase
                        VECCOPY(verts[i].v, id->V[i]);
                }
        }
+       
        return 1;
 }