Initial commit of cloth modifier from branch rev 13453
authorDaniel Genrich <daniel.genrich@gmx.net>
Tue, 29 Jan 2008 21:01:12 +0000 (21:01 +0000)
committerDaniel Genrich <daniel.genrich@gmx.net>
Tue, 29 Jan 2008 21:01:12 +0000 (21:01 +0000)
36 files changed:
CMakeLists.txt
SConstruct
extern/SConscript
extern/bullet2/src/Bullet-C-Api.h [new file with mode: 0644]
extern/bullet2/src/BulletDynamics/CMakeLists.txt
extern/bullet2/src/BulletDynamics/Dynamics/Bullet-C-API.cpp [new file with mode: 0644]
extern/bullet2/src/SConscript
intern/elbeem/CMakeLists.txt
intern/elbeem/SConscript
intern/elbeem/intern/isosurface.cpp
source/Makefile
source/blender/blenkernel/BKE_cloth.h [new file with mode: 0644]
source/blender/blenkernel/BKE_collision.h [new file with mode: 0644]
source/blender/blenkernel/BKE_modifier.h
source/blender/blenkernel/CMakeLists.txt
source/blender/blenkernel/SConscript
source/blender/blenkernel/intern/Makefile
source/blender/blenkernel/intern/cloth.c [new file with mode: 0644]
source/blender/blenkernel/intern/collision.c [new file with mode: 0644]
source/blender/blenkernel/intern/implicit.c [new file with mode: 0644]
source/blender/blenkernel/intern/kdop.c [new file with mode: 0644]
source/blender/blenkernel/intern/modifier.c
source/blender/blenloader/intern/readfile.c
source/blender/blenloader/intern/writefile.c
source/blender/include/butspace.h
source/blender/makesdna/DNA_cloth_types.h [new file with mode: 0644]
source/blender/makesdna/DNA_modifier_types.h
source/blender/makesdna/intern/makesdna.c
source/blender/src/buttons_editing.c
source/blender/src/buttons_object.c
source/blender/src/editmesh.c
source/blender/src/editobject.c
source/blender/src/transform_conversions.c
source/blender/src/transform_generics.c
source/blender/src/vpaint.c
tools/btools.py

index 1804eb5444663163c025fb7e4675387e49b7a19a..905fb6a26739080304d3d697bbc9e37d282b5ef8 100644 (file)
@@ -65,8 +65,9 @@ OPTION(WITH_ELBEEM            "Enable Elbeem (Fluid Simulation)"                      ON)
 OPTION(WITH_QUICKTIME          "Enable Quicktime Support"                              OFF)
 OPTION(WITH_OPENEXR            "Enable OpenEXR Support (http://www.openexr.com)"       OFF)
 OPTION(WITH_FFMPEG             "Enable FFMPeg Support (http://ffmpeg.mplayerhq.hu/)"   OFF)
-OPTION(WITH_OPENAL             "Enable OpenAL Support (http://www.openal.org)" ON)
-OPTION(YESIAMSTUPID            "Enable execution on 64-bit platforms"                                  OFF)
+OPTION(WITH_OPENAL             "Enable OpenAL Support (http://www.openal.org)"         ON)
+OPTION(YESIAMSTUPID            "Enable execution on 64-bit platforms"                  OFF)
+OPTION(WITH_OPENMP             "Enable OpenMP (has to be supported by the compiler)"   OFF)
 
 IF(NOT WITH_GAMEENGINE AND WITH_PLAYER)
   MESSAGE("WARNING: WITH_PLAYER needs WITH_GAMEENGINE")
@@ -184,6 +185,13 @@ IF(UNIX)
 
   SET(LLIBS "-lXi -lutil -lc -lm -lpthread -lstdc++")
 
+  IF(WITH_OPENMP)
+    SET(LLIBS "${LLIBS} -lgomp ")
+    SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp ")
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp ")
+  ENDIF(WITH_OPENMP)
+
+
   SET(PLATFORM_CFLAGS "-pipe -fPIC -funsigned-char -fno-strict-aliasing -DXP_UNIX -Wno-char-subscripts")
 
   SET(PLATFORM_LINKFLAGS "-pthread")
@@ -270,6 +278,11 @@ IF(WIN32)
   SET(CMAKE_C_FLAGS_MINSIZEREL "/D_CRT_NONSTDC_NO_DEPRECATE /D_CRT_SECURE_NO_DEPRECATE /D_SCL_SECURE_NO_DEPRECATE /wd4800 /wd4244 /wd4305 /O1 /Ob1 /DNDEBUG /EHsc /MT /W3 /nologo /J" CACHE STRING "MSVC MT flags " FORCE)
   SET(CMAKE_C_FLAGS_RELWITHDEBINFO "/D_CRT_NONSTDC_NO_DEPRECATE /D_CRT_SECURE_NO_DEPRECATE /D_SCL_SECURE_NO_DEPRECATE /wd4800 /wd4244 /wd4305 /O2 /Ob1 /DNDEBUG /EHsc /MT /W3 /nologo /Zi /J" CACHE STRING "MSVC MT flags " FORCE)
 
+  IF(WITH_OPENMP)
+    SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} /openmp ")
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /openmp ")
+  ENDIF(WITH_OPENMP)
+
   SET(SDL ${LIBDIR}/sdl)
   SET(SDL_INC ${SDL}/include)
   SET(SDL_LIB SDL)
@@ -347,6 +360,12 @@ IF(APPLE)
   SET(PLATFORM_CFLAGS "-pipe -fPIC -funsigned-char -fno-strict-aliasing")
   SET(PLATFORM_LINKFLAGS "-fexceptions -framework CoreServices -framework Foundation -framework IOKit -framework AppKit -framework Carbon -framework AGL -framework AudioUnit -framework AudioToolbox -framework CoreAudio -framework QuickTime")
 
+  IF(WITH_OPENMP)
+    SET(LLIBS "${LLIBS} -lgomp ")
+    SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp ")
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp ")
+  ENDIF(WITH_OPENMP)
+
   SET(SDL ${LIBDIR}/sdl)
   SET(SDL_INC ${SDL}/include)
   SET(SDL_LIB SDL)
index 51f925f634b3598f15b0c8dac8bc6e526b41e2bf..1886b6486b8ebc90eb5804bf05813fb8cefe5596 100644 (file)
@@ -179,6 +179,18 @@ if env['BF_NO_ELBEEM'] == 1:
     env['CXXFLAGS'].append('-DDISABLE_ELBEEM')
     env['CCFLAGS'].append('-DDISABLE_ELBEEM')
 
+if env['WITH_BF_OPENMP'] == 1:
+       if env['OURPLATFORM']=='win32-vc':
+               env.Append(LINKFLAGS=['/openmp'])
+               env['CCFLAGS'].append('/openmp')
+               env['CPPFLAGS'].append('/openmp')
+               env['CXXFLAGS'].append('/openmp')
+       else:
+               env.Append(LINKFLAGS=['-lgomp'])
+               env['CCFLAGS'].append('-fopenmp')
+               env['CPPFLAGS'].append('-fopenmp')
+               env['CXXFLAGS'].append('-fopenmp')
+
 #check for additional debug libnames
 
 if env.has_key('BF_DEBUG_LIBS'):
index eefd1bbb635f2b5c591245220a052d6cf01f9054..bcf67055788c7aca7a72cad5c1e1b4049761034d 100644 (file)
@@ -5,8 +5,9 @@ Import('env')
 if env['WITH_BF_GAMEENGINE']:
     SConscript(['qhull/SConscript',
             'solid/SConscript'])
-    if env['WITH_BF_BULLET']:
-        SConscript(['bullet2/src/SConscript'])
+
+if env['WITH_BF_BULLET']:
+    SConscript(['bullet2/src/SConscript'])
 
 if env['WITH_BF_INTERNATIONAL']:
     SConscript(['bFTGL/SConscript'])
diff --git a/extern/bullet2/src/Bullet-C-Api.h b/extern/bullet2/src/Bullet-C-Api.h
new file mode 100644 (file)
index 0000000..078dcae
--- /dev/null
@@ -0,0 +1,37 @@
+/*
+Bullet Continuous Collision Detection and Physics Library
+Copyright (c) 2003-2006 Erwin Coumans  http://continuousphysics.com/Bullet/
+
+This software is provided 'as-is', without any express or implied warranty.
+In no event will the authors be held liable for any damages arising from the use of this software.
+Permission is granted to anyone to use this software for any purpose, 
+including commercial applications, and to alter it and redistribute it freely, 
+subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required.
+2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+3. This notice may not be removed or altered from any source distribution.
+*/
+
+/*
+       Draft high-level generic physics C-API. For low-level access, use the physics SDK native API's.
+       Work in progress, functionality will be added on demand.
+
+       If possible, use the richer Bullet C++ API, by including "btBulletDynamicsCommon.h"
+*/
+
+#ifndef BULLET_C_API_H
+#define BULLET_C_API_H
+
+#ifdef __cplusplus
+extern "C" { 
+#endif
+
+double plNearestPoints(float p1[3], float p2[3], float p3[3], float q1[3], float q2[3], float q3[3], float *pa, float *pb, float normal[3]);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif //BULLET_C_API_H
+
index 79e07b7f77bca7ab6664bcc26f6e4a583cc813ae..8598575799a64733ff580d9f15d90929360d92ef 100644 (file)
@@ -13,6 +13,7 @@ ADD_LIBRARY(LibBulletDynamics
        ConstraintSolver/btTypedConstraint.cpp
        Dynamics/btDiscreteDynamicsWorld.cpp
        Dynamics/btSimpleDynamicsWorld.cpp
+       Dynamics/Bullet-C-API.cpp
        Dynamics/btRigidBody.cpp
        Vehicle/btRaycastVehicle.cpp
        Vehicle/btWheelInfo.cpp
diff --git a/extern/bullet2/src/BulletDynamics/Dynamics/Bullet-C-API.cpp b/extern/bullet2/src/BulletDynamics/Dynamics/Bullet-C-API.cpp
new file mode 100644 (file)
index 0000000..d4f3761
--- /dev/null
@@ -0,0 +1,115 @@
+/*
+Bullet Continuous Collision Detection and Physics Library
+Copyright (c) 2003-2006 Erwin Coumans  http://continuousphysics.com/Bullet/
+
+This software is provided 'as-is', without any express or implied warranty.
+In no event will the authors be held liable for any damages arising from the use of this software.
+Permission is granted to anyone to use this software for any purpose, 
+including commercial applications, and to alter it and redistribute it freely, 
+subject to the following restrictions:
+
+1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required.
+2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software.
+3. This notice may not be removed or altered from any source distribution.
+*/
+
+/*
+       Draft high-level generic physics C-API. For low-level access, use the physics SDK native API's.
+       Work in progress, functionality will be added on demand.
+
+       If possible, use the richer Bullet C++ API, by including <src/btBulletDynamicsCommon.h>
+*/
+
+#include "Bullet-C-Api.h"
+#include "btBulletDynamicsCommon.h"
+#include "LinearMath/btAlignedAllocator.h"
+
+
+#include "LinearMath/btVector3.h"
+#include "LinearMath/btScalar.h"       
+#include "LinearMath/btMatrix3x3.h"
+#include "LinearMath/btTransform.h"
+#include "BulletCollision/NarrowPhaseCollision/btVoronoiSimplexSolver.h"
+#include "BulletCollision/CollisionShapes/btTriangleShape.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
+#include "BulletCollision/NarrowPhaseCollision/btPointCollector.h"
+#include "BulletCollision/NarrowPhaseCollision/btVoronoiSimplexSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btSubSimplexConvexCast.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpaPenetrationDepthSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpa.h"
+#include "BulletCollision/CollisionShapes/btMinkowskiSumShape.h"
+#include "BulletCollision/NarrowPhaseCollision/btSubSimplexConvexCast.h"
+#include "BulletCollision/NarrowPhaseCollision/btVoronoiSimplexSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
+
+#include "BulletCollision/NarrowPhaseCollision/btDiscreteCollisionDetectorInterface.h"
+#include "BulletCollision/NarrowPhaseCollision/btSimplexSolverInterface.h"
+#include "BulletCollision/NarrowPhaseCollision/btMinkowskiPenetrationDepthSolver.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkPairDetector.h"
+#include "BulletCollision/NarrowPhaseCollision/btGjkEpaPenetrationDepthSolver.h"
+#include "LinearMath/btStackAlloc.h"
+
+extern "C"
+double plNearestPoints(float p1[3], float p2[3], float p3[3], float q1[3], float q2[3], float q3[3], float *pa, float *pb, float normal[3])
+{
+       btTriangleShape trishapeA(btVector3(p1[0], p1[1], p1[2]), btVector3(p2[0], p2[1], p2[2]), btVector3(p3[0], p3[1], p3[2]));
+       trishapeA.setMargin(0.000001f);
+       
+       btTriangleShape trishapeB(btVector3(q1[0], q1[1], q1[2]), btVector3(q2[0], q2[1], q2[2]), btVector3(q3[0], q3[1], q3[2]));
+       trishapeB.setMargin(0.000001f);
+       
+       // btVoronoiSimplexSolver sGjkSimplexSolver;
+       // btGjkEpaPenetrationDepthSolver penSolverPtr; 
+       
+       static btSimplexSolverInterface sGjkSimplexSolver;
+       sGjkSimplexSolver.reset();
+       
+       static btGjkEpaPenetrationDepthSolver Solver0;
+       static btMinkowskiPenetrationDepthSolver Solver1;
+               
+       btConvexPenetrationDepthSolver* Solver = NULL;
+       
+       Solver = &Solver1;      
+               
+       btGjkPairDetector convexConvex(&trishapeA ,&trishapeB,&sGjkSimplexSolver,Solver);
+       
+       convexConvex.m_catchDegeneracies = 1;
+       
+       // btGjkPairDetector convexConvex(&trishapeA ,&trishapeB,&sGjkSimplexSolver,0);
+       
+       btPointCollector gjkOutput;
+       btGjkPairDetector::ClosestPointInput input;
+       
+       btStackAlloc gStackAlloc(1024*1024*2);
+       input.m_stackAlloc = &gStackAlloc;
+       
+       btTransform tr;
+       tr.setIdentity();
+       
+       input.m_transformA = tr;
+       input.m_transformB = tr;
+       
+       convexConvex.getClosestPoints(input, gjkOutput, 0);
+       
+       
+       if (gjkOutput.m_hasResult)
+       {
+               
+               pb[0] = pa[0] = gjkOutput.m_pointInWorld[0];
+               pb[1] = pa[1] = gjkOutput.m_pointInWorld[1];
+               pb[2] = pa[2] = gjkOutput.m_pointInWorld[2];
+
+               pb[0]+= gjkOutput.m_normalOnBInWorld[0] * gjkOutput.m_distance;
+               pb[1]+= gjkOutput.m_normalOnBInWorld[1] * gjkOutput.m_distance;
+               pb[2]+= gjkOutput.m_normalOnBInWorld[2] * gjkOutput.m_distance;
+               
+               normal[0] = gjkOutput.m_normalOnBInWorld[0];
+               normal[1] = gjkOutput.m_normalOnBInWorld[1];
+               normal[2] = gjkOutput.m_normalOnBInWorld[2];
+
+               return gjkOutput.m_distance;
+       }
+       return -1.0f;   
+}
index 6280c49066d1987777df352d5ae21d1686ba5980..19702782b0dd5bd55a85d16bd419558d9c572b71 100644 (file)
@@ -34,6 +34,7 @@ bulletdyn_src = ["BulletDynamics/ConstraintSolver/btContactConstraint.cpp",
                  "BulletDynamics/Dynamics/btSimpleDynamicsWorld.cpp",
                  "BulletDynamics/Dynamics/btRigidBody.cpp",
                  "BulletDynamics/Vehicle/btRaycastVehicle.cpp",
+                "BulletDynamics/Dynamics/Bullet-C-API.cpp",
                  "BulletDynamics/Vehicle/btWheelInfo.cpp"]
 collision_src = ["BulletCollision/BroadphaseCollision/btAxisSweep3.cpp",
                  "BulletCollision/BroadphaseCollision/btBroadphaseProxy.cpp",
index 86a60180b053225eb81c1625687e8649be07878d..ac35b8ce00f04c2d483afecd6c38f9972eedd731 100644 (file)
@@ -36,5 +36,9 @@ IF(WINDOWS)
     ADD_DEFINITIONS(-DUSE_MSVC6FIXES)
 ENDIF(WINDOWS)
 
+IF(WITH_OPENMP)
+    ADD_DEFINITIONS(-DPARALLEL)
+ENDIF(WITH_OPENMP)
+
 BLENDERLIB_NOLIST(bf_elbeem "${SRC}" "${INC}")
 #, libtype='blender', priority=0 )
index bb6637ba32d1109cf2cd431cb6edbe3eba1b5a85..bdcb050798731480fef9681b3f77f408d57303d6 100644 (file)
@@ -5,7 +5,11 @@ Import('env')
 
 sources = env.Glob('intern/*.cpp')
 
-defs = 'NOGUI ELBEEM_BLENDER=1'
+defs = ' NOGUI ELBEEM_BLENDER=1'
+
+if env['WITH_BF_OPENMP'] == 1:
+    defs += ' PARALLEL'
+
 if env['OURPLATFORM']=='win32-vc':
     defs += ' USE_MSVC6FIXES'
 incs = env['BF_PNG_INC'] + ' ' + env['BF_ZLIB_INC'] + ' ' +env['BF_SDL_INC']
index 1b0ba13c70710b55d4874bbb61de44ed0d542ba4..f90621f3b73cc6871a7c929e4f6522dead052a3b 100644 (file)
 #define round(x) (x)
 #endif
 
+#if PARALLEL==1        
+#include <omp.h>
+#endif
+
 /******************************************************************************
  * Constructor
  *****************************************************************************/
@@ -160,13 +164,6 @@ void IsoSurface::triangulate( void )
                mpEdgeVerticesZ[i] = -1;
        }
 
-       ntlVec3Gfx pos[8];
-       float value[8];
-       int cubeIndex;      // index entry of the cube 
-       int triIndices[12]; // vertex indices 
-       int *eVert[12];
-       IsoLevelVertex ilv;
-
        // edges between which points?
        const int mcEdges[24] = { 
                0,1,  1,2,  2,3,  3,0,
@@ -193,7 +190,12 @@ void IsoSurface::triangulate( void )
                                px = mStart[0]-gsx*0.5;
                                for(int i=1;i<(mSizex-2);i++) {
                                        px += gsx;
-
+                                       int cubeIndex;      // index entry of the cube 
+                                       float value[8];
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       
                                        value[0] = *getData(i  ,j  ,k  );
                                        value[1] = *getData(i+1,j  ,k  );
                                        value[2] = *getData(i+1,j+1,k  );
@@ -239,6 +241,7 @@ void IsoSurface::triangulate( void )
                                        eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
 
                                        // grid positions
+                                       ntlVec3Gfx pos[8];
                                        pos[0] = ntlVec3Gfx(px    ,py    ,pz);
                                        pos[1] = ntlVec3Gfx(px+gsx,py    ,pz);
                                        pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz);
@@ -344,10 +347,7 @@ void IsoSurface::triangulate( void )
                if(mUseFullEdgeArrays) {
                        errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!");
                }
-
-               // subdiv local arrays
-               gfxReal orgval[8];
-               gfxReal subdAr[2][11][11]; // max 10 subdivs!
+               
                ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex];
 
                // construct pointers
@@ -408,13 +408,25 @@ void IsoSurface::triangulate( void )
 
                debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
                pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
+       
                for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
                        pz += gsz;
                        const int k = ok/mSubdivs;
                        if(k<=0) continue; // skip zero plane
+#if PARALLEL==1        
+#pragma omp parallel for
+#endif 
                        for(int j=1;j<(mSizey-2);j++) {
                                for(int i=1;i<(mSizex-2);i++) {
-
+                                       float value[8];
+                                       ntlVec3Gfx pos[8];
+                                       int cubeIndex;      // index entry of the cube 
+                                       int triIndices[12]; // vertex indices 
+                                       int *eVert[12];
+                                       IsoLevelVertex ilv;
+                                       gfxReal orgval[8];
+                                       gfxReal subdAr[2][11][11]; // max 10 subdivs!
+                                       
                                        orgval[0] = *getData(i  ,j  ,k  );
                                        orgval[1] = *getData(i+1,j  ,k  );
                                        orgval[2] = *getData(i+1,j+1,k  ); // with subdivs
@@ -426,6 +438,7 @@ void IsoSurface::triangulate( void )
 
                                        // prebuild subsampled array slice
                                        const int sdkOffset = ok-k*mSubdivs; 
+
                                        for(int sdk=0; sdk<2; sdk++) 
                                                for(int sdj=0; sdj<mSubdivs+1; sdj++) 
                                                        for(int sdi=0; sdi<mSubdivs+1; sdi++) {
@@ -580,8 +593,13 @@ void IsoSurface::triangulate( void )
 
                                                                                // init isolevel vertex
                                                                                ilv.v = p1 + (p2-p1)*mu; // with subdivs
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                                               {
                                                                                mPoints.push_back( ilv );
                                                                                triIndices[e] = (mPoints.size()-1);
+                                                                               }
                                                                                // store vertex 
                                                                                *eVert[ e ] = triIndices[e]; 
                                                                        }       else {
@@ -591,23 +609,27 @@ void IsoSurface::triangulate( void )
                                                                } // along all edges 
                                                        }
                                                        // removed cutoff treatment...
-
+                                                       
                                                        // Create the triangles... 
+#if PARALLEL==1        
+#pragma omp critical
+#endif
+                                                       {
                                                        for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] );
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); // with subdivs
                                                                mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
                                                                //errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
-                                                       }
-
                                                        } // triangles in edge table?
+                                                       }
+                                                       }
                                                        
                                                }//si
                                        }// sj
 
                                }//i
                        }// j
-
+                       
                        // copy edge arrays
                        for(int j=0;j<(mSizey-0)*mSubdivs;j++) 
                                for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
index 3553c723b5179ea67ec7cfc5f67e4960a7a1403b..5c72940d9a2fa91656f21d82ff442040210b7eae 100644 (file)
@@ -151,9 +151,11 @@ ifneq ($(NAN_NO_KETSJI),true)
     COMLIB += $(OCGDIR)/gameengine/ketsji/KXNetwork/$(DEBUG_DIR)libKXNetwork.a
     COMLIB += $(OCGDIR)/gameengine/Network/$(DEBUG_DIR)libNetwork.a
     COMLIB += $(OCGDIR)/gameengine/Network/LoopBackNetwork/$(DEBUG_DIR)libLoopBackNetwork.a
-    COMLIB += $(NAN_BULLET2)/lib/libbullet2.a
 endif
 
+# Required by cloth, not gameengine only anymore
+COMLIB += $(NAN_BULLET2)/lib/$(DEBUG_DIR)libbullet2.a
+
 COMLIB += $(NAN_GUARDEDALLOC)/lib/libguardedalloc.a
 COMLIB += $(NAN_MEMUTIL)/lib/libmemutil.a
 COMLIB += $(NAN_BMFONT)/lib/$(DEBUG_DIR)libbmfont.a
diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h
new file mode 100644 (file)
index 0000000..530601f
--- /dev/null
@@ -0,0 +1,256 @@
+/**
+ * BKE_cloth.h
+ *
+ * $Id: BKE_cloth.h,v 1.1 2007/08/01 02:07:27 daniel Exp $
+ *
+ * ***** BEGIN GPL/BL DUAL 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.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+#ifndef BKE_CLOTH_H
+#define BKE_CLOTH_H
+
+#include "float.h"
+#include "BLI_linklist.h"
+#include "BKE_collision.h"
+#include "BKE_customdata.h"
+#include "BKE_DerivedMesh.h"
+#include "DNA_cloth_types.h"
+#include "DNA_customdata_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_object_types.h"
+
+struct Object;
+struct Cloth;
+struct MFace;
+struct DerivedMesh;
+struct ClothModifierData;
+struct CollisionTree;
+
+// this is needed for inlining behaviour
+#ifndef _WIN32
+#define LINUX
+#define DO_INLINE inline
+#else
+#define DO_INLINE
+#endif
+
+#define CLOTH_MAX_THREAD 2
+
+/**
+ * The definition of a cloth vertex.
+ */
+typedef struct ClothVertex
+{
+       int     flags;          /* General flags per vertex.            */
+       float   v [3];          /* The velocity of the point.           */
+       float   xconst [3];     /* constrained position                 */
+       float   x [3];          /* The current position of this vertex. */
+       float   xold [3];       /* The previous position of this vertex.*/
+       float   tx [3];         /* temporary position */
+       float   txold [3];      /* temporary old position */
+       float   tv[3];          /* temporary "velocity", mostly used as tv = tx-txold */
+       float   mass;           /* mass / weight of the vertex          */
+       float   goal;           /* goal, from SB                        */
+       float   impulse[3];     /* used in collision.c */
+       unsigned int impulse_count; /* same as above */
+       float   avg_spring_len; /* average length of connected springs, UNUSED ATM */
+       float   struct_stiff;
+       float   bend_stiff;
+       float   shear_stiff;
+}
+ClothVertex;
+
+/**
+ * The definition of a spring.
+ */
+typedef struct ClothSpring
+{
+       int     ij;             /* Pij from the paper, one end of the spring.   */
+       int     kl;             /* Pkl from the paper, one end of the spring.   */
+       float   restlen;        /* The original length of the spring.   */
+       int     matrix_index;   /* needed for implicit solver (fast lookup) */
+       int     type;           /* types defined in BKE_cloth.h ("springType") */
+       int     flags;          /* defined in BKE_cloth.h, e.g. deactivated due to tearing */
+       float dfdx[3][3];
+       float dfdv[3][3];
+       float f[3];
+       float   stiffness;      /* stiffness factor from the vertex groups */
+}
+ClothSpring;
+
+/* goal defines */
+#define SOFTGOALSNAP  0.999f
+
+/* This is approximately the smallest number that can be
+* represented by a float, given its precision. */
+#define ALMOST_ZERO            FLT_EPSILON
+
+// some macro enhancements for vector treatment
+#define VECADDADD(v1,v2,v3)    {*(v1)+= *(v2) + *(v3); *(v1+1)+= *(v2+1) + *(v3+1); *(v1+2)+= *(v2+2) + *(v3+2);}
+#define VECSUBADD(v1,v2,v3)    {*(v1)-= *(v2) + *(v3); *(v1+1)-= *(v2+1) + *(v3+1); *(v1+2)-= *(v2+2) + *(v3+2);}
+#define VECADDSUB(v1,v2,v3)    {*(v1)+= *(v2) - *(v3); *(v1+1)+= *(v2+1) - *(v3+1); *(v1+2)+= *(v2+2) - *(v3+2);}
+#define VECSUBADDSS(v1,v2,aS,v3,bS)    {*(v1)-= *(v2)*aS + *(v3)*bS; *(v1+1)-= *(v2+1)*aS + *(v3+1)*bS; *(v1+2)-= *(v2+2)*aS + *(v3+2)*bS;}
+#define VECADDSUBSS(v1,v2,aS,v3,bS)    {*(v1)+= *(v2)*aS - *(v3)*bS; *(v1+1)+= *(v2+1)*aS - *(v3+1)*bS; *(v1+2)+= *(v2+2)*aS - *(v3+2)*bS;}
+#define VECADDSS(v1,v2,aS,v3,bS)       {*(v1)= *(v2)*aS + *(v3)*bS; *(v1+1)= *(v2+1)*aS + *(v3+1)*bS; *(v1+2)= *(v2+2)*aS + *(v3+2)*bS;}
+#define VECADDS(v1,v2,v3,bS)   {*(v1)= *(v2) + *(v3)*bS; *(v1+1)= *(v2+1) + *(v3+1)*bS; *(v1+2)= *(v2+2) + *(v3+2)*bS;}
+#define VECSUBMUL(v1,v2,aS)    {*(v1)-= *(v2) * aS; *(v1+1)-= *(v2+1) * aS; *(v1+2)-= *(v2+2) * aS;}
+#define VECSUBS(v1,v2,v3,bS)   {*(v1)= *(v2) - *(v3)*bS; *(v1+1)= *(v2+1) - *(v3+1)*bS; *(v1+2)= *(v2+2) - *(v3+2)*bS;}
+#define VECSUBSB(v1,v2, v3,bS)         {*(v1)= (*(v2)- *(v3))*bS; *(v1+1)= (*(v2+1) - *(v3+1))*bS; *(v1+2)= (*(v2+2) - *(v3+2))*bS;}
+#define VECMULS(v1,aS)         {*(v1)*= aS; *(v1+1)*= aS; *(v1+2)*= *aS;}
+#define VECADDMUL(v1,v2,aS)    {*(v1)+= *(v2) * aS; *(v1+1)+= *(v2+1) * aS; *(v1+2)+= *(v2+2) * aS;}
+
+/* SIMULATION FLAGS: goal flags,.. */
+/* These are the bits used in SimSettings.flags. */
+typedef enum
+{
+       CLOTH_SIMSETTINGS_FLAG_RESET = ( 1 << 1 ),      // The CM object requires a reinitializaiton.
+                                        CLOTH_SIMSETTINGS_FLAG_COLLOBJ = ( 1 << 2 ),// object is only collision object, no cloth simulation is done
+                                                        CLOTH_SIMSETTINGS_FLAG_GOAL = ( 1 << 3 ),      // we have goals enabled
+                                                                        CLOTH_SIMSETTINGS_FLAG_TEARING = ( 1 << 4 ),// true if tearing is enabled
+                                                                               CLOTH_SIMSETTINGS_FLAG_CCACHE_PROTECT = ( 1 << 5 ), // true if tearing is enabled
+                                                                               CLOTH_SIMSETTINGS_FLAG_EDITMODE = ( 1 << 6 ), // are we in editmode? -several things disabled
+                                                                               CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE = (1 << 7), /* force cache freeing */
+                                                                               CLOTH_SIMSETTINGS_FLAG_SCALING = (1 << 8), /* is advanced scaling active? */
+                                                                               CLOTH_SIMSETTINGS_FLAG_LOADED = (1 << 9), /* did we just got load? */
+} CLOTH_SIMSETTINGS_FLAGS;
+
+/* COLLISION FLAGS */
+typedef enum
+{
+       CLOTH_COLLSETTINGS_FLAG_ENABLED = ( 1 << 1 ), /* enables cloth - object collisions */
+                                           CLOTH_COLLSETTINGS_FLAG_SELF = ( 1 << 2 ), /* unused */
+} CLOTH_COLLISIONSETTINGS_FLAGS;
+
+/* Spring types as defined in the paper.*/
+typedef enum
+{
+       CLOTH_SPRING_TYPE_STRUCTURAL = 0,
+ CLOTH_SPRING_TYPE_SHEAR,
+ CLOTH_SPRING_TYPE_BENDING,
+} CLOTH_SPRING_TYPES;
+
+/* SPRING FLAGS */
+typedef enum
+{
+       CLOTH_SPRING_FLAG_DEACTIVATE = ( 1 << 1 ),
+                                        CLOTH_SPRING_FLAG_NEEDED = ( 1 << 2 ), // springs has values to be applied
+} CLOTH_SPRINGS_FLAGS;
+
+/* Bits to or into the ClothVertex.flags. */
+#define CLOTH_VERT_FLAG_PINNED 1
+#define CLOTH_VERT_FLAG_COLLISION 2
+
+typedef void ( *CM_COLLISION_RESPONSE ) ( ClothModifierData *clmd, CollisionModifierData *collmd, CollisionTree *tree1, CollisionTree *tree2 );
+
+
+/////////////////////////////////////////////////
+// collision.c
+////////////////////////////////////////////////
+
+// needed for implicit.c
+void bvh_collision_response ( ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree * tree1, CollisionTree * tree2 );
+int cloth_bvh_objcollision ( ClothModifierData * clmd, float step, float dt );
+
+int bvh_traverse ( ClothModifierData * clmd, CollisionModifierData * collmd, CollisionTree * tree1, CollisionTree * tree2, float step, CM_COLLISION_RESPONSE collision_response );
+////////////////////////////////////////////////
+
+
+////////////////////////////////////////////////
+// implicit.c
+////////////////////////////////////////////////
+               
+// needed for cloth.c
+int implicit_init ( Object *ob, ClothModifierData *clmd );
+int implicit_free ( ClothModifierData *clmd );
+int implicit_solver ( Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors );
+void implicit_set_positions ( ClothModifierData *clmd );
+////////////////////////////////////////////////
+
+
+/////////////////////////////////////////////////
+// cloth.c
+////////////////////////////////////////////////
+
+// needed for modifier.c
+void cloth_free_modifier_extern (ClothModifierData *clmd);
+void cloth_free_modifier (Object *ob, ClothModifierData *clmd);
+void cloth_init (ClothModifierData *clmd);
+DerivedMesh *clothModifier_do(ClothModifierData *clmd,Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc);
+
+void cloth_update_normals (ClothVertex *verts, int nVerts, MFace *face, int totface);
+
+// needed for collision.c
+void bvh_update_from_cloth(ClothModifierData *clmd, int moving);
+
+// needed for editmesh.c
+void cloth_write_cache(Object *ob, ClothModifierData *clmd, float framenr);
+int cloth_read_cache(Object *ob, ClothModifierData *clmd, float framenr);
+
+// needed for button_object.c
+void cloth_clear_cache(Object *ob, ClothModifierData *clmd, float framenr);
+
+////////////////////////////////////////////////
+
+
+/* Typedefs for function pointers we need for solvers and collision detection. */
+typedef void ( *CM_COLLISION_SELF ) ( ClothModifierData *clmd, int step );
+typedef void ( *CM_COLLISION_OBJ ) ( ClothModifierData *clmd, int step, CM_COLLISION_RESPONSE collision_response );
+
+
+/* This enum provides the IDs for our solvers. */
+// only one available in the moment
+typedef enum {
+       CM_IMPLICIT = 0,
+} CM_SOLVER_ID;
+
+
+/* This structure defines how to call the solver.
+*/
+typedef struct
+{
+       char            *name;
+       CM_SOLVER_ID    id;
+       int     ( *init ) ( Object *ob, ClothModifierData *clmd );
+       int     ( *solver ) ( Object *ob, float framenr, ClothModifierData *clmd, ListBase *effectors );
+       int     ( *free ) ( ClothModifierData *clmd );
+}
+CM_SOLVER_DEF;
+
+/* used for caching in implicit.c */
+typedef struct Frame
+{
+       ClothVertex *verts;
+       ClothSpring *springs;
+       unsigned int numverts, numsprings;
+       float time; /* we need float since we want to support sub-frames */
+}
+Frame;
+
+#endif
+
diff --git a/source/blender/blenkernel/BKE_collision.h b/source/blender/blenkernel/BKE_collision.h
new file mode 100644 (file)
index 0000000..aa46bc6
--- /dev/null
@@ -0,0 +1,176 @@
+/**
+ * BKE_cloth.h
+ *
+ * $Id: BKE_cloth.h,v 1.1 2007/08/01 02:07:27 daniel Exp $
+ *
+ * ***** BEGIN GPL/BL DUAL 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.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+#ifndef BKE_COLLISIONS_H
+#define BKE_COLLISIONS_H
+
+#include <math.h>
+#include "float.h"
+#include <stdlib.h>
+#include <string.h>
+
+/* types */
+#include "BLI_linklist.h"
+#include "BKE_collision.h"
+#include "BKE_customdata.h"
+#include "BKE_DerivedMesh.h"
+#include "DNA_cloth_types.h"
+#include "DNA_customdata_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_object_types.h"
+
+struct Object;
+struct Cloth;
+struct MFace;
+struct DerivedMesh;
+struct ClothModifierData;
+struct CollisionTree;
+
+
+////////////////////////////////////////
+// used in kdop.c and collision.c
+////////////////////////////////////////
+typedef struct CollisionTree
+{
+       struct CollisionTree *nodes[4]; // 4 children --> quad-tree
+       struct CollisionTree *parent;
+       struct CollisionTree *nextLeaf;
+       struct CollisionTree *prevLeaf;
+       float   bv[26]; // Bounding volume of all nodes / we have 7 axes on a 14-DOP
+       unsigned int tri_index; // this saves the index of the face
+       // int point_index[4]; // supports up to 4 points in a leaf
+       int     count_nodes; // how many nodes are used
+       int     traversed;  // how many nodes already traversed until this level?
+       int     isleaf;
+}
+CollisionTree;
+
+typedef struct BVH
+{
+       unsigned int    numfaces;
+       unsigned int    numverts;
+       MVert           *current_x; // e.g. txold in clothvertex
+       MVert           *current_xold; // e.g. tx in clothvertex
+       MFace           *mfaces; // just a pointer to the original datastructure
+       struct LinkNode *tree;
+       CollisionTree   *root; // TODO: saving the root --> is this really needed? YES!
+       CollisionTree   *leaf_tree; /* Tail of the leaf linked list.    */
+       CollisionTree   *leaf_root;     /* Head of the leaf linked list.        */
+       float           epsilon; /* epslion is used for inflation of the k-dop     */
+       int             flags; /* bvhFlags */
+}
+BVH;
+////////////////////////////////////////
+
+
+////////////////////////////////////////
+// used for collisions in kdop.c and also collision.c
+////////////////////////////////////////
+/* used for collisions in collision.c */
+typedef struct CollPair
+{
+       unsigned int face1; // cloth face
+       unsigned int face2; // object face
+       double distance; // magnitude of vector
+       float normal[3];
+       float vector[3]; // unnormalized collision vector: p2-p1
+       float pa[3], pb[3]; // collision point p1 on face1, p2 on face2
+       int lastsign; // indicates if the distance sign has changed, unused itm
+       float time; // collision time, from 0 up to 1
+       unsigned int ap1, ap2, ap3, bp1, bp2, bp3, bp4;
+       unsigned int pointsb[4];
+}
+CollPair;
+
+/* used for collisions in collision.c */
+typedef struct EdgeCollPair
+{
+       unsigned int p11, p12, p21, p22;
+       float normal[3];
+       float vector[3];
+       float time;
+       int lastsign;
+       float pa[3], pb[3]; // collision point p1 on face1, p2 on face2
+}
+EdgeCollPair;
+
+/* used for collisions in collision.c */
+typedef struct FaceCollPair
+{
+       unsigned int p11, p12, p13, p21;
+       float normal[3];
+       float vector[3];
+       float time;
+       int lastsign;
+       float pa[3], pb[3]; // collision point p1 on face1, p2 on face2
+}
+FaceCollPair;
+////////////////////////////////////////
+
+
+
+/////////////////////////////////////////////////
+// forward declarations
+/////////////////////////////////////////////////
+
+// NOTICE: mvert-routines for building + update the BVH are the most native ones
+
+// builds bounding volume hierarchy
+void bvh_build (BVH *bvh);
+BVH *bvh_build_from_mvert (MFace *mfaces, unsigned int numfaces, MVert *x, unsigned int numverts, float epsilon);
+
+// frees the same
+void bvh_free ( BVH * bvh );
+
+// checks two bounding volume hierarchies for potential collisions and returns some list with those
+
+
+// update bounding volumes, needs updated positions in bvh->x
+void bvh_update_from_mvert(BVH * bvh, MVert *x, unsigned int numverts, MVert *xnew, int moving);
+void bvh_update(BVH * bvh, int moving);
+
+LinkNode *BLI_linklist_append_fast ( LinkNode **listp, void *ptr );
+
+// move Collision modifier object inter-frame with step = [0,1]
+// defined in collisions.c
+void collision_move_object(CollisionModifierData *collmd, float step, float prevstep);
+
+// interface for collision functions
+void collisions_compute_barycentric (float pv[3], float p1[3], float p2[3], float p3[3], float *w1, float *w2, float *w3);
+void interpolateOnTriangle(float to[3], float v1[3], float v2[3], float v3[3], double w1, double w2, double w3);
+
+/////////////////////////////////////////////////
+
+#endif
+
index 4b679c1587f658acde977b2a8b7348004c93d978..fba30264feadf87573cfdd97e98fa6da75e9b6bc 100644 (file)
@@ -287,7 +287,9 @@ int           modifiers_getCageIndex(struct Object *ob,
                                      int *lastPossibleCageIndex_r);
 
 int           modifiers_isSoftbodyEnabled(struct Object *ob);
+int           modifiers_isClothEnabled(struct Object *ob);
 int           modifiers_isParticleEnabled(struct Object *ob);
+
 struct Object *modifiers_isDeformedByArmature(struct Object *ob);
 struct Object *modifiers_isDeformedByLattice(struct Object *ob);
 int           modifiers_usesArmature(struct Object *ob, struct bArmature *arm);
index 0b87f0c1d987a80ef8334dde3c6c27d5cc474218..30f21ef83cc6a4acae1b747b53d616e6bd685fd0 100644 (file)
@@ -34,7 +34,7 @@ SET(INC
   ../python ../render/extern/include ../../../intern/decimation/extern
   ../imbuf ../avi ../../../intern/elbeem/extern ../../../intern/opennl/extern
   ../../../intern/iksolver/extern ../blenloader ../quicktime
-  ../../../intern/bmfont
+  ../../../intern/bmfont ../../../extern/bullet2/src
   ../nodes
   ${SDL_INC}
   ${ZLIB_INC}
index 744ccf2749197efce49b37312fb32f60439348dc..f8f2f0b9f57f89be09faa7bd60b63ce905c39ee9 100644 (file)
@@ -7,6 +7,7 @@ incs = '. #/intern/guardedalloc ../include ../blenlib ../makesdna'
 incs += ' ../python ../render/extern/include #/intern/decimation/extern'
 incs += ' ../imbuf ../avi #/intern/elbeem/extern ../nodes'
 incs += ' #/intern/iksolver/extern ../blenloader ../quicktime'
+incs += ' #/extern/bullet2/src'
 incs += ' #/intern/bmfont'
 incs += ' #/intern/opennl/extern'
 
index 3d37f11a051d2d767e4003dc85a274c6a7b85909..488d8801e418bb480ac3a8444ba3994c2c598295 100644 (file)
@@ -80,6 +80,9 @@ CPPFLAGS += -I../../nodes
 # path to our own external headerfiles
 CPPFLAGS += -I..
 
+# path to bullet2, for cloth
+CPPFLAGS += -I../../../../extern/bullet2/src
+
 ifeq ($(WITH_FREETYPE2), true)
     CPPFLAGS += -DWITH_FREETYPE2
     CPPFLAGS += -I$(NAN_FREETYPE)/include
diff --git a/source/blender/blenkernel/intern/cloth.c b/source/blender/blenkernel/intern/cloth.c
new file mode 100644 (file)
index 0000000..f871a30
--- /dev/null
@@ -0,0 +1,1441 @@
+/*  cloth.c
+*
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) Blender Foundation
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): none yet.
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "MEM_guardedalloc.h"
+
+/* types */
+#include "DNA_curve_types.h"
+#include "DNA_object_types.h"
+#include "DNA_object_force.h"
+#include "DNA_cloth_types.h"
+#include "DNA_key_types.h"
+#include "DNA_mesh_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_edgehash.h"
+#include "BLI_linklist.h"
+
+#include "BKE_curve.h"
+#include "BKE_deform.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_displist.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+#include "BKE_key.h"
+#include "BKE_mesh.h"
+#include "BKE_object.h"
+#include "BKE_cloth.h"
+#include "BKE_modifier.h"
+#include "BKE_utildefines.h"
+#include "BKE_DerivedMesh.h"
+#include "BIF_editdeform.h"
+#include "BIF_editkey.h"
+#include "DNA_screen_types.h"
+#include "BSE_headerbuttons.h"
+#include "BIF_screen.h"
+#include "BIF_space.h"
+#include "mydevice.h"
+
+#include "BKE_pointcache.h"
+
+#ifdef _WIN32
+void tstart ( void )
+{}
+void tend ( void )
+{
+}
+double tval()
+{
+       return 0;
+}
+#else
+#include <sys/time.h>
+                        static struct timeval _tstart, _tend;
+        static struct timezone tz;
+        void tstart ( void )
+{
+       gettimeofday ( &_tstart, &tz );
+}
+void tend ( void )
+{
+       gettimeofday ( &_tend,&tz );
+}
+double tval()
+{
+       double t1, t2;
+       t1 = ( double ) _tstart.tv_sec + ( double ) _tstart.tv_usec/ ( 1000*1000 );
+       t2 = ( double ) _tend.tv_sec + ( double ) _tend.tv_usec/ ( 1000*1000 );
+       return t2-t1;
+}
+#endif
+
+/* Our available solvers. */
+// 255 is the magic reserved number, so NEVER try to put 255 solvers in here!
+// 254 = MAX!
+static CM_SOLVER_DEF   solvers [] =
+{
+       { "Implicit", CM_IMPLICIT, implicit_init, implicit_solver, implicit_free },
+        // { "Implicit C++", CM_IMPLICITCPP, implicitcpp_init, implicitcpp_solver, implicitcpp_free },
+};
+
+/* ********** cloth engine ******* */
+/* Prototypes for internal functions.
+*/
+static void cloth_to_object (Object *ob,  ClothModifierData *clmd, DerivedMesh *dm);
+static void cloth_from_mesh ( Object *ob, ClothModifierData *clmd, DerivedMesh *dm );
+static int cloth_from_object(Object *ob, ClothModifierData *clmd, DerivedMesh *dm, float framenr);
+int cloth_build_springs ( ClothModifierData *clmd, DerivedMesh *dm );
+static void cloth_apply_vgroup ( ClothModifierData *clmd, DerivedMesh *dm );
+
+
+/******************************************************************************
+*
+* External interface called by modifier.c clothModifier functions.
+*
+******************************************************************************/
+/**
+ * cloth_init -  creates a new cloth simulation.
+ *
+ * 1. create object
+ * 2. fill object with standard values or with the GUI settings if given
+ */
+void cloth_init ( ClothModifierData *clmd )
+{      
+       /* Initialize our new data structure to reasonable values. */
+       clmd->sim_parms->gravity [0] = 0.0;
+       clmd->sim_parms->gravity [1] = 0.0;
+       clmd->sim_parms->gravity [2] = -9.81;
+       clmd->sim_parms->structural = 30.0;
+       clmd->sim_parms->shear = 30.0;
+       clmd->sim_parms->bending = 1.0;
+       clmd->sim_parms->Cdis = 5.0;
+       clmd->sim_parms->Cvi = 1.0;
+       clmd->sim_parms->mass = 1.0f;
+       clmd->sim_parms->stepsPerFrame = 5;
+       clmd->sim_parms->sim_time = 1.0;
+       clmd->sim_parms->flags = 0;
+       clmd->sim_parms->solver_type = 0;
+       clmd->sim_parms->preroll = 0;
+       clmd->sim_parms->maxspringlen = 10;
+       clmd->sim_parms->firstframe = 1;
+       clmd->sim_parms->lastframe = 250;
+       clmd->sim_parms->vgroup_mass = 0;
+       clmd->sim_parms->lastcachedframe = 0;
+       clmd->sim_parms->editedframe = 0;
+       clmd->sim_parms->autoprotect = 25;
+       clmd->sim_parms->firstcachedframe = -1.0;
+       
+       clmd->coll_parms->self_friction = 5.0;
+       clmd->coll_parms->friction = 10.0;
+       clmd->coll_parms->loop_count = 3;
+       clmd->coll_parms->epsilon = 0.015f;
+       clmd->coll_parms->flags = CLOTH_COLLSETTINGS_FLAG_ENABLED;
+
+       /* These defaults are copied from softbody.c's
+       * softbody_calc_forces() function.
+       */
+       clmd->sim_parms->eff_force_scale = 1000.0;
+       clmd->sim_parms->eff_wind_scale = 250.0;
+
+       // also from softbodies
+       clmd->sim_parms->maxgoal = 1.0f;
+       clmd->sim_parms->mingoal = 0.0f;
+       clmd->sim_parms->defgoal = 0.0f;
+       clmd->sim_parms->goalspring = 100.0f;
+       clmd->sim_parms->goalfrict = 0.0f;
+}
+
+
+BVH *bvh_build_from_cloth (ClothModifierData *clmd, float epsilon)
+{
+       unsigned int i = 0;
+       BVH     *bvh=NULL;
+       Cloth *cloth = clmd->clothObject;
+       ClothVertex *verts = NULL;
+
+       if(!clmd)
+               return NULL;
+
+       cloth = clmd->clothObject;
+
+       if(!cloth)
+               return NULL;
+       
+       verts = cloth->verts;
+       
+       bvh = MEM_callocN(sizeof(BVH), "BVH");
+       if (bvh == NULL) 
+       {
+               printf("bvh: Out of memory.\n");
+               return NULL;
+       }
+       
+       // springs = cloth->springs;
+       // numsprings = cloth->numsprings;
+       
+       bvh->flags = 0;
+       bvh->leaf_tree = NULL;
+       bvh->leaf_root = NULL;
+       bvh->tree = NULL;
+
+       bvh->epsilon = epsilon;
+       bvh->numfaces = cloth->numfaces;
+       bvh->mfaces = cloth->mfaces;
+
+       bvh->numverts = cloth->numverts;
+       
+       bvh->current_x = MEM_callocN ( sizeof ( MVert ) * bvh->numverts, "bvh->current_x" );
+       bvh->current_xold = MEM_callocN ( sizeof ( MVert ) * bvh->numverts, "bvh->current_xold" );
+       
+       for(i = 0; i < bvh->numverts; i++)
+       {
+               VECCOPY(bvh->current_x[i].co, verts[i].tx);
+               VECCOPY(bvh->current_xold[i].co, verts[i].txold);
+       }
+       
+       bvh_build (bvh);
+       
+       return bvh;
+}
+
+void bvh_update_from_cloth(ClothModifierData *clmd, int moving)
+{
+       unsigned int i = 0;
+       Cloth *cloth = clmd->clothObject;
+       BVH *bvh = cloth->tree;
+       ClothVertex *verts = cloth->verts;
+       
+       if(!bvh)
+               return;
+       
+       if(cloth->numverts!=bvh->numverts)
+               return;
+       
+       if(cloth->verts)
+       {
+               for(i = 0; i < bvh->numverts; i++)
+               {
+                       VECCOPY(bvh->current_x[i].co, verts[i].tx);
+                       VECCOPY(bvh->current_xold[i].co, verts[i].txold);
+               }
+       }
+       
+       bvh_update(bvh, moving);
+}
+
+// unused in the moment, cloth needs quads from mesh
+DerivedMesh *CDDM_convert_to_triangle ( DerivedMesh *dm )
+{
+       DerivedMesh *result = NULL;
+       int i;
+       int numverts = dm->getNumVerts ( dm );
+       int numedges = dm->getNumEdges ( dm );
+       int numfaces = dm->getNumFaces ( dm );
+
+       MVert *mvert = CDDM_get_verts ( dm );
+       MEdge *medge = CDDM_get_edges ( dm );
+       MFace *mface = CDDM_get_faces ( dm );
+
+       MVert *mvert2;
+       MFace *mface2;
+       unsigned int numtris=0;
+       unsigned int numquads=0;
+       int a = 0;
+       int random = 0;
+       int firsttime = 0;
+       float vec1[3], vec2[3], vec3[3], vec4[3], vec5[3];
+       float mag1=0, mag2=0;
+
+       for ( i = 0; i < numfaces; i++ )
+       {
+               if ( mface[i].v4 )
+                       numquads++;
+               else
+                       numtris++;
+       }
+
+       result = CDDM_from_template ( dm, numverts, 0, numtris + 2*numquads );
+
+       if ( !result )
+               return NULL;
+
+       // do verts
+       mvert2 = CDDM_get_verts ( result );
+       for ( a=0; a<numverts; a++ )
+       {
+               MVert *inMV;
+               MVert *mv = &mvert2[a];
+
+               inMV = &mvert[a];
+
+               DM_copy_vert_data ( dm, result, a, a, 1 );
+               *mv = *inMV;
+       }
+
+
+       // do faces
+       mface2 = CDDM_get_faces ( result );
+       for ( a=0, i=0; a<numfaces; a++ )
+       {
+               MFace *mf = &mface2[i];
+               MFace *inMF;
+               inMF = &mface[a];
+
+               /*
+               DM_copy_face_data(dm, result, a, i, 1);
+
+               *mf = *inMF;
+               */
+
+               if ( mface[a].v4 && random==1 )
+               {
+                       mf->v1 = mface[a].v2;
+                       mf->v2 = mface[a].v3;
+                       mf->v3 = mface[a].v4;
+               }
+               else
+               {
+                       mf->v1 = mface[a].v1;
+                       mf->v2 = mface[a].v2;
+                       mf->v3 = mface[a].v3;
+               }
+
+               mf->v4 = 0;
+               mf->flag |= ME_SMOOTH;
+
+               test_index_face ( mf, NULL, 0, 3 );
+
+               if ( mface[a].v4 )
+               {
+                       MFace *mf2;
+
+                       i++;
+
+                       mf2 = &mface2[i];
+                       /*
+                       DM_copy_face_data(dm, result, a, i, 1);
+
+                       *mf2 = *inMF;
+                       */
+
+                       if ( random==1 )
+                       {
+                               mf2->v1 = mface[a].v1;
+                               mf2->v2 = mface[a].v2;
+                               mf2->v3 = mface[a].v4;
+                       }
+                       else
+                       {
+                               mf2->v1 = mface[a].v4;
+                               mf2->v2 = mface[a].v1;
+                               mf2->v3 = mface[a].v3;
+                       }
+                       mf2->v4 = 0;
+                       mf2->flag |= ME_SMOOTH;
+
+                       test_index_face ( mf2, NULL, 0, 3 );
+               }
+
+               i++;
+       }
+
+       CDDM_calc_edges ( result );
+       CDDM_calc_normals ( result );
+
+       return result;
+
+}
+
+
+DerivedMesh *CDDM_create_tearing ( ClothModifierData *clmd, DerivedMesh *dm )
+{
+       DerivedMesh *result = NULL;
+       unsigned int i = 0, a = 0, j=0;
+       int numverts = dm->getNumVerts ( dm );
+       int numedges = dm->getNumEdges ( dm );
+       int numfaces = dm->getNumFaces ( dm );
+
+       MVert *mvert = CDDM_get_verts ( dm );
+       MEdge *medge = CDDM_get_edges ( dm );
+       MFace *mface = CDDM_get_faces ( dm );
+
+       MVert *mvert2;
+       MFace *mface2;
+       unsigned int numtris=0;
+       unsigned int numquads=0;
+       EdgeHash *edgehash = NULL;
+       Cloth *cloth = clmd->clothObject;
+       ClothSpring *springs = cloth->springs;
+       unsigned int numsprings = cloth->numsprings;
+
+       // create spring tearing hash
+       edgehash = BLI_edgehash_new();
+
+       for ( i = 0; i < numsprings; i++ )
+       {
+               if ( ( springs[i].flags & CLOTH_SPRING_FLAG_DEACTIVATE )
+                                    && ( !BLI_edgehash_haskey ( edgehash, springs[i].ij, springs[i].kl ) ) )
+               {
+                       BLI_edgehash_insert ( edgehash, springs[i].ij, springs[i].kl, NULL );
+                       BLI_edgehash_insert ( edgehash, springs[i].kl, springs[i].ij, NULL );
+                       j++;
+               }
+       }
+
+       // printf("found %d tears\n", j);
+
+       result = CDDM_from_template ( dm, numverts, 0, numfaces );
+
+       if ( !result )
+               return NULL;
+
+       // do verts
+       mvert2 = CDDM_get_verts ( result );
+       for ( a=0; a<numverts; a++ )
+       {
+               MVert *inMV;
+               MVert *mv = &mvert2[a];
+
+               inMV = &mvert[a];
+
+               DM_copy_vert_data ( dm, result, a, a, 1 );
+               *mv = *inMV;
+       }
+
+
+       // do faces
+       mface2 = CDDM_get_faces ( result );
+       for ( a=0, i=0; a<numfaces; a++ )
+       {
+               MFace *mf = &mface2[i];
+               MFace *inMF;
+               inMF = &mface[a];
+
+               /*
+               DM_copy_face_data(dm, result, a, i, 1);
+
+               *mf = *inMF;
+               */
+
+               if ( ( !BLI_edgehash_haskey ( edgehash, mface[a].v1, mface[a].v2 ) )
+                                     && ( !BLI_edgehash_haskey ( edgehash, mface[a].v2, mface[a].v3 ) )
+                                     && ( !BLI_edgehash_haskey ( edgehash, mface[a].v3, mface[a].v4 ) )
+                                     && ( !BLI_edgehash_haskey ( edgehash, mface[a].v4, mface[a].v1 ) ) )
+               {
+                       mf->v1 = mface[a].v1;
+                       mf->v2 = mface[a].v2;
+                       mf->v3 = mface[a].v3;
+                       mf->v4 = mface[a].v4;
+
+                       test_index_face ( mf, NULL, 0, 4 );
+
+                       i++;
+               }
+       }
+
+       CDDM_lower_num_faces ( result, i );
+       CDDM_calc_edges ( result );
+       CDDM_calc_normals ( result );
+
+       BLI_edgehash_free ( edgehash, NULL );
+
+       return result;
+}
+
+int modifiers_indexInObject(Object *ob, ModifierData *md_seek);
+
+int cloth_read_cache(Object *ob, ClothModifierData *clmd, float framenr)
+{
+       FILE *fp = NULL;
+       int stack_index = -1;
+       unsigned int a, ret = 1;
+       Cloth *cloth = clmd->clothObject;
+       
+       if(!cloth)
+               return 0;
+       
+       stack_index = modifiers_indexInObject(ob, (ModifierData *)clmd);
+       
+       fp = BKE_ptcache_id_fopen((ID *)ob, 'r', framenr, stack_index);
+       if(!fp)
+               ret = 0;
+       else {
+               for(a = 0; a < cloth->numverts; a++)
+               {
+                       if(fread(&cloth->verts[a].x, sizeof(float), 3, fp) != 3) 
+                       {
+                               ret = 0;
+                               break;
+                       }
+                       if(fread(&cloth->verts[a].xconst, sizeof(float), 3, fp) != 3) 
+                       {
+                               ret = 0;
+                               break;
+                       }
+                       if(fread(&cloth->verts[a].v, sizeof(float), 3, fp) != 3) 
+                       {
+                               ret = 0;
+                               break;
+                       }
+               }
+               
+               fclose(fp);
+               
+               if(clmd->sim_parms->lastcachedframe < framenr)
+               {
+                       if(G.rt > 0)
+                               printf("cloth_read_cache problem: lnex - f#: %f, lastCF: %d\n", framenr, clmd->sim_parms->lastcachedframe);
+               }
+       }
+       
+       if(G.rt > 0)
+               printf("cloth_read_cache: %f\n", framenr);
+       
+       return ret;
+}
+
+void cloth_clear_cache(Object *ob, ClothModifierData *clmd, float framenr)
+{
+       int stack_index = -1;
+       
+       // don't do anything as long as we're in editmode!
+       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_EDITMODE)
+       {
+               /* delete cache free request */
+               clmd->sim_parms->flags &= ~CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE;
+               
+               return;
+       }
+       
+       /* clear cache if specific frame cleaning requested or cache is not protected */
+       if((!(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_CCACHE_PROTECT)) || (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE))
+       {
+               stack_index = modifiers_indexInObject(ob, (ModifierData *)clmd);
+               
+               BKE_ptcache_id_clear((ID *)ob, PTCACHE_CLEAR_AFTER, framenr, stack_index);
+               
+               /* update last cached frame # */
+               clmd->sim_parms->lastcachedframe = framenr;
+               
+               /* update first cached frame # */
+               if((framenr < clmd->sim_parms->firstcachedframe) && (clmd->sim_parms->firstcachedframe >=0.0))
+                       clmd->sim_parms->firstcachedframe = -1.0;
+               
+               if(G.rt > 0)
+                       printf("cloth_clear_cache: %f\n", framenr);
+       }
+       
+       /* delete cache free request */
+       clmd->sim_parms->flags &= ~CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE;
+       
+       
+}
+void cloth_write_cache(Object *ob, ClothModifierData *clmd, float framenr)
+{
+       FILE *fp = NULL;
+       int stack_index = -1;
+       unsigned int a;
+       Cloth *cloth = clmd->clothObject;
+       
+       if(G.rt > 0)
+               printf("cloth_write_cache: %f\n", framenr);
+       
+       if(!cloth)
+       {
+               if(G.rt > 0)
+                       printf("cloth_write_cache: no cloth\n");
+               return;
+       }
+       
+       stack_index = modifiers_indexInObject(ob, (ModifierData *)clmd);
+       
+       fp = BKE_ptcache_id_fopen((ID *)ob, 'w', framenr, stack_index);
+       if(!fp)
+       {
+               if(G.rt > 0)
+                       printf("cloth_write_cache: no fp\n");
+               return;
+       }
+       
+       for(a = 0; a < cloth->numverts; a++)
+       {
+               fwrite(&cloth->verts[a].x, sizeof(float),3,fp);
+               fwrite(&cloth->verts[a].xconst, sizeof(float),3,fp);
+               fwrite(&cloth->verts[a].v, sizeof(float),3,fp);
+       }
+       
+       /* update last cached frame # */
+       clmd->sim_parms->lastcachedframe = MAX2(clmd->sim_parms->lastcachedframe, framenr);
+       
+       /* update first cached frame # */
+       if((clmd->sim_parms->firstcachedframe < 0.0) || ((framenr < clmd->sim_parms->firstcachedframe) && (clmd->sim_parms->firstcachedframe > 0.0)))
+               clmd->sim_parms->firstcachedframe = framenr;
+       
+       if(G.rt > 0)
+               printf("lcf: %d, framenr: %f\n", clmd->sim_parms->lastcachedframe, framenr);
+
+       fclose(fp);
+}
+
+
+
+/************************************************
+ * clothModifier_do - main simulation function
+************************************************/
+DerivedMesh *clothModifier_do(ClothModifierData *clmd,Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
+
+{
+       unsigned int i;
+       Cloth *cloth = clmd->clothObject;
+       float framenr = G.scene->r.cfra;
+       float current_time = bsystem_time ( ob, ( float ) G.scene->r.cfra, 0.0 );
+       ListBase *effectors = NULL;
+       ClothVertex *verts = NULL;
+       float deltaTime = current_time - clmd->sim_parms->sim_time;
+       unsigned int numverts = -1;
+       unsigned int numedges = -1;
+       unsigned int numfaces = -1;
+       MVert *mvert = NULL;
+       MEdge *medge = NULL;
+       MFace *mface = NULL;
+       DerivedMesh *result = NULL;
+       
+       if(G.rt > 0)
+               printf("clothModifier_do start\n");
+       
+       /* we're getting called two times during file load,
+       resulting in a not valid G.relbase on the first time (cache makes problems)
+       --> just return back */
+       if((clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_LOADED)&& (!G.relbase_valid))
+       {
+               clmd->sim_parms->flags &= ~CLOTH_SIMSETTINGS_FLAG_LOADED;
+               return dm;
+       }
+       
+       result = CDDM_copy(dm);
+       
+       if(!result)
+       {
+               return dm;
+       }
+       
+       numverts = result->getNumVerts(result);
+       numedges = result->getNumEdges(result);
+       numfaces = result->getNumFaces(result);
+       mvert = dm->getVertArray(result);
+       medge = dm->getEdgeArray(result);
+       mface = dm->getFaceArray(result);
+       
+       /* check if cache is active / if file is already saved */
+       /*
+       if ((!G.relbase_valid) && ( deltaTime != 1.0f ))
+       {
+       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESET;
+}
+       */
+       
+       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_RESET)
+       {       
+               cloth_free_modifier (ob, clmd);
+               if(G.rt > 0)
+                       printf("clothModifier_do CLOTH_SIMSETTINGS_FLAG_RESET\n");
+       }
+       
+       // unused in the moment, calculated seperately in implicit.c
+       clmd->sim_parms->dt = 1.0f / clmd->sim_parms->stepsPerFrame;
+       
+       if ( ( clmd->clothObject == NULL ) || (clmd->clothObject && (numverts != clmd->clothObject->numverts )) )
+       {
+               /* only force free the cache if we have a different number of verts */
+               if(clmd->clothObject && (numverts != clmd->clothObject->numverts ))
+               {
+                       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE;
+                       cloth_free_modifier ( ob, clmd );
+               }
+               
+               cloth_clear_cache(ob, clmd, 0);
+                               
+               if ( !cloth_from_object ( ob, clmd, result, framenr ) )
+                       return result;
+       
+               if ( clmd->clothObject == NULL )
+                       return result;
+       
+               cloth = clmd->clothObject;
+               
+               if(!cloth_read_cache(ob, clmd, framenr))
+               {
+                       /* save first frame in case we have a reseted object 
+                       and we move one frame forward.
+                       In that case we would only start with the SECOND frame
+                       if we don't save the current state before 
+                       TODO PROBLEM: IMHO we can't track external movement from the
+                       first frame in this case! */
+                       /*
+                       if ( deltaTime == 1.0f )
+                       cloth_write_cache(ob, clmd, framenr-1.0);
+                       */
+                       if(G.rt > 0)
+                               printf("cloth_from_object NO cloth_read_cache cloth_write_cache\n");
+               }
+               else
+               {
+                       if(G.rt > 0)
+                               printf("cloth_from_object cloth_read_cache\n");
+                       
+                       implicit_set_positions(clmd);
+               }
+               
+               clmd->sim_parms->sim_time = current_time;
+       }
+       
+       // only be active during a specific period:
+       // that's "first frame" and "last frame" on GUI
+       
+       // TODO: enable later again after refactoring
+       if ( current_time < clmd->sim_parms->firstframe )
+       {
+               return result;
+       }
+       else if ( current_time > clmd->sim_parms->lastframe )
+       {
+               int stack_index = modifiers_indexInObject(ob, (ModifierData *)clmd);
+                       
+               if(BKE_ptcache_id_exist((ID *)ob, clmd->sim_parms->lastcachedframe, stack_index))
+               {
+                       if(cloth_read_cache(ob, clmd, clmd->sim_parms->lastcachedframe))
+                       {
+                               implicit_set_positions(clmd);
+                               
+                               // Copy the result back to the object.
+                               cloth_to_object (ob, clmd, result);
+                       }
+               }
+               return result;
+       }
+       
+       /* nice moving one frame forward */
+       if ( deltaTime == 1.0f )
+       {
+               clmd->sim_parms->sim_time = current_time;
+               
+               if(G.rt > 0)
+                       printf("clothModifier_do deltaTime=1\n");
+               
+               if(!cloth_read_cache(ob, clmd, framenr))
+               {
+                       verts = cloth->verts;
+
+                       // Force any pinned verts to their constrained location.
+                       for ( i = 0; i < clmd->clothObject->numverts; i++, verts++ )
+                       {
+                               // Save the previous position.
+                               VECCOPY ( verts->xold, verts->xconst );
+                               VECCOPY ( verts->txold, verts->x );
+
+                               // Get the current position.
+                               VECCOPY ( verts->xconst, mvert[i].co );
+                               Mat4MulVecfl ( ob->obmat, verts->xconst );
+                       }
+                       
+                       tstart();
+
+                       // Call the solver.
+                       if ( solvers [clmd->sim_parms->solver_type].solver )
+                               solvers [clmd->sim_parms->solver_type].solver ( ob, framenr, clmd, effectors );
+
+                       tend();
+                       // printf ( "Cloth simulation time: %f\n", ( float ) tval() );
+
+                       cloth_write_cache(ob, clmd, framenr);
+                       
+                       // check for autoprotection
+                       if(framenr >= clmd->sim_parms->autoprotect)
+                       {
+                               if(G.rt > 0)
+                                       printf("fr#: %f, auto: %d\n", framenr, clmd->sim_parms->autoprotect);
+                               clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_CCACHE_PROTECT;
+                       }
+                       if(G.rt > 0)
+                               printf("clothModifier_do deltaTime=1 cachewrite\n");
+               }
+               else
+               {
+                       if(G.rt > 0)
+                               printf("clothModifier_do deltaTime=1 cacheread\n");
+                       implicit_set_positions(clmd);
+               }
+               
+               // Copy the result back to the object.
+               cloth_to_object (ob, clmd, result);
+       }
+       else if(deltaTime == 0.0f) 
+       {       
+               if(G.rt > 0)
+                       printf("clothModifier_do deltaTime!=1 clmd->clothObject != NULL\n");
+               if(cloth_read_cache(ob, clmd, framenr))
+               {
+                       cloth_to_object (ob, clmd, result);
+                       implicit_set_positions(clmd);
+               }
+               else /* same cache parts are missing */
+               {
+                       /*
+                       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESET;
+                       */
+                       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_CCACHE_FFREE;
+                       cloth_clear_cache(ob, clmd, 0);
+                       
+                       cloth_write_cache(ob, clmd, framenr);
+               }
+       }
+       else
+       {       
+               if(G.rt > 0)
+                       printf("clothModifier_do deltaTime!=1 clmd->clothObject != NULL\n");
+               if(cloth_read_cache(ob, clmd, framenr))
+               {
+                       cloth_to_object (ob, clmd, result);
+                       implicit_set_positions(clmd);
+                       clmd->sim_parms->sim_time = current_time;
+               }
+               else
+               {
+                       /* jump to a non-existing frame makes sim reset */
+                       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESET;
+               }
+       }
+       
+       return result;
+}
+
+/* frees all */
+void cloth_free_modifier ( Object *ob, ClothModifierData *clmd )
+{
+       Cloth   *cloth = NULL;
+       
+       if ( !clmd )
+               return;
+
+       cloth = clmd->clothObject;
+
+       
+       if ( cloth )
+       {       
+               // If our solver provides a free function, call it
+               if ( solvers [clmd->sim_parms->solver_type].free )
+               {
+                       solvers [clmd->sim_parms->solver_type].free ( clmd );
+               }
+
+               // Free the verts.
+               if ( cloth->verts != NULL )
+                       MEM_freeN ( cloth->verts );
+
+               cloth->verts = NULL;
+               cloth->numverts = 0;
+
+               // Free the springs.
+               if ( cloth->springs != NULL )
+               {
+                       LinkNode *search = cloth->springs;
+                       while(search)
+                       {
+                               ClothSpring *spring = search->link;
+                                               
+                               MEM_freeN ( spring );
+                               search = search->next;
+                       }
+                       BLI_linklist_free(cloth->springs, NULL);
+               
+                       cloth->springs = NULL;
+               }
+
+               cloth->springs = NULL;
+               cloth->numsprings = 0;
+
+               // free BVH collision tree
+               if ( cloth->tree )
+                       bvh_free ( ( BVH * ) cloth->tree );
+
+               // we save our faces for collision objects
+               if ( cloth->mfaces )
+                       MEM_freeN ( cloth->mfaces );
+               /*
+               if(clmd->clothObject->facemarks)
+               MEM_freeN(clmd->clothObject->facemarks);
+               */
+               MEM_freeN ( cloth );
+               clmd->clothObject = NULL;
+       }
+       clmd->sim_parms->flags &= ~CLOTH_SIMSETTINGS_FLAG_RESET;
+}
+
+/* frees all */
+void cloth_free_modifier_extern ( ClothModifierData *clmd )
+{
+       Cloth   *cloth = NULL;
+       if(G.rt > 0)
+               printf("cloth_free_modifier_extern\n");
+       
+       if ( !clmd )
+               return;
+
+       cloth = clmd->clothObject;
+       
+       if ( cloth )
+       {       
+               if(G.rt > 0)
+                       printf("cloth_free_modifier_extern in\n");
+               
+               // If our solver provides a free function, call it
+               if ( solvers [clmd->sim_parms->solver_type].free )
+               {
+                       solvers [clmd->sim_parms->solver_type].free ( clmd );
+               }
+
+               // Free the verts.
+               if ( cloth->verts != NULL )
+                       MEM_freeN ( cloth->verts );
+
+               cloth->verts = NULL;
+               cloth->numverts = 0;
+
+               // Free the springs.
+               if ( cloth->springs != NULL )
+               {
+                       LinkNode *search = cloth->springs;
+                       while(search)
+                       {
+                               ClothSpring *spring = search->link;
+                                               
+                               MEM_freeN ( spring );
+                               search = search->next;
+                       }
+                       BLI_linklist_free(cloth->springs, NULL);
+               
+                       cloth->springs = NULL;
+               }
+
+               cloth->springs = NULL;
+               cloth->numsprings = 0;
+
+               // free BVH collision tree
+               if ( cloth->tree )
+                       bvh_free ( ( BVH * ) cloth->tree );
+
+               // we save our faces for collision objects
+               if ( cloth->mfaces )
+                       MEM_freeN ( cloth->mfaces );
+               /*
+               if(clmd->clothObject->facemarks)
+               MEM_freeN(clmd->clothObject->facemarks);
+               */
+               MEM_freeN ( cloth );
+               clmd->clothObject = NULL;
+       }
+}
+
+/******************************************************************************
+*
+* Internal functions.
+*
+******************************************************************************/
+
+/**
+ * cloth_to_object - copies the deformed vertices to the object.
+ *
+ **/
+static void cloth_to_object (Object *ob,  ClothModifierData *clmd, DerivedMesh *dm)
+{
+       unsigned int    i = 0;
+       MVert *mvert = NULL;
+       unsigned int numverts;
+       Cloth *cloth = clmd->clothObject;
+
+       if (clmd->clothObject) {
+               /* inverse matrix is not uptodate... */
+               Mat4Invert (ob->imat, ob->obmat);
+
+               mvert = CDDM_get_verts(dm);
+               numverts = dm->getNumVerts(dm);
+
+               for (i = 0; i < numverts; i++)
+               {
+                       VECCOPY (mvert[i].co, cloth->verts[i].x);
+                       Mat4MulVecfl (ob->imat, mvert[i].co);   /* cloth is in global coords */
+               }
+       }
+}
+
+
+/**
+ * cloth_apply_vgroup - applies a vertex group as specified by type
+ *
+ **/
+/* can be optimized to do all groups in one loop */
+static void cloth_apply_vgroup ( ClothModifierData *clmd, DerivedMesh *dm )
+{
+       unsigned int i = 0;
+       unsigned int j = 0;
+       MDeformVert *dvert = NULL;
+       Cloth *clothObj = NULL;
+       unsigned int numverts = dm->getNumVerts ( dm );
+       float goalfac = 0;
+       ClothVertex *verts = NULL;
+       // clmd->sim_parms->vgroup_mass
+
+       clothObj = clmd->clothObject;
+
+       if ( !dm )
+               return;
+       
+       numverts = dm->getNumVerts ( dm );
+
+       verts = clothObj->verts;
+       
+       if (((clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_SCALING ) || 
+                    (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL )) && 
+                    ((clmd->sim_parms->vgroup_mass>0) || 
+                    (clmd->sim_parms->vgroup_struct>0)||
+                    (clmd->sim_parms->vgroup_bend>0)))
+       {
+               for ( i = 0; i < numverts; i++, verts++ )
+               {
+                       dvert = dm->getVertData ( dm, i, CD_MDEFORMVERT );
+                       if ( dvert )
+                       {
+                               for ( j = 0; j < dvert->totweight; j++ )
+                               {
+                                       if (( dvert->dw[j].def_nr == (clmd->sim_parms->vgroup_mass-1)) && (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL ))
+                                       {
+                                               verts->goal = dvert->dw [j].weight;
+                                               goalfac= 1.0f;
+                                               
+                                               /*
+                                               // Kicking goal factor to simplify things...who uses that anyway?
+                                               // ABS ( clmd->sim_parms->maxgoal - clmd->sim_parms->mingoal );
+                                               */
+                                               
+                                               verts->goal  = ( float ) pow ( verts->goal , 4.0f );
+                                               if ( verts->goal >=SOFTGOALSNAP )
+                                                       verts->flags |= CLOTH_VERT_FLAG_PINNED;
+                                       }
+                                       
+                                       if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_SCALING )
+                                       {
+                                               if( dvert->dw[j].def_nr == (clmd->sim_parms->vgroup_struct-1))
+                                               {
+                                                       verts->struct_stiff = dvert->dw [j].weight;
+                                                       verts->shear_stiff = dvert->dw [j].weight;
+                                               }
+                                               
+                                               if( dvert->dw[j].def_nr == (clmd->sim_parms->vgroup_bend-1))
+                                               {
+                                                       verts->bend_stiff = dvert->dw [j].weight;
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+}
+
+/*
+helper function to get proper spring length
+when object is rescaled
+*/
+float cloth_globallen ( float *v1,float *v2,Object *ob )
+{
+       float p1[3],p2[3];
+       VECCOPY ( p1,v1 );
+       Mat4MulVecfl ( ob->obmat, p1 );
+       VECCOPY ( p2,v2 );
+       Mat4MulVecfl ( ob->obmat, p2 );
+       return VecLenf ( p1,p2 );
+}
+
+static int cloth_from_object(Object *ob, ClothModifierData *clmd, DerivedMesh *dm, float framenr)
+{
+       unsigned int i = 0;
+       MVert *mvert = NULL;
+       ClothVertex *verts = NULL;
+       float tnull[3] = {0,0,0};
+       int cache_there = 0;
+
+       // If we have a clothObject, free it. 
+       if ( clmd->clothObject != NULL )
+       {
+               cloth_free_modifier ( ob, clmd );
+               if(G.rt > 0)
+                       printf("cloth_free_modifier cloth_from_object\n");
+       }
+
+       // Allocate a new cloth object.
+       clmd->clothObject = MEM_callocN ( sizeof ( Cloth ), "cloth" );
+       if ( clmd->clothObject )
+       {
+               clmd->clothObject->old_solver_type = 255;
+               // clmd->clothObject->old_collision_type = 255;
+       }
+       else if ( !clmd->clothObject )
+       {
+               modifier_setError ( & ( clmd->modifier ), "Out of memory on allocating clmd->clothObject." );
+               return 0;
+       }
+
+       // mesh input objects need DerivedMesh
+       if ( !dm )
+               return 0;
+
+       cloth_from_mesh ( ob, clmd, dm );
+
+       if((clmd->sim_parms->firstcachedframe < 0.0) || ((clmd->sim_parms->firstcachedframe >= 0.0) && (!cloth_read_cache(ob, clmd, clmd->sim_parms->firstcachedframe))))
+       {
+               // no cache there
+               cache_there = 0;
+               if(G.rt > 0)
+                       printf("cache_there = 0\n");
+       }
+       else
+       {
+               // we have a cache
+               cache_there = 1;
+               if(G.rt > 0)
+                       printf("cache_there = 1, fcf: %d\n", clmd->sim_parms->firstcachedframe);
+       }
+       
+       // create springs 
+       clmd->clothObject->springs = NULL;
+       clmd->clothObject->numsprings = -1;
+       
+       mvert = dm->getVertArray ( dm );
+       verts = clmd->clothObject->verts;
+
+       // set initial values
+       for ( i = 0; i < dm->getNumVerts(dm); i++, verts++ )
+       {
+               if(!cache_there)
+               {
+                       VECCOPY ( verts->x, mvert[i].co );
+                       Mat4MulVecfl ( ob->obmat, verts->x );
+               }
+               
+               verts->mass = clmd->sim_parms->mass;
+
+               if ( clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL )
+                       verts->goal= clmd->sim_parms->defgoal;
+               else
+                       verts->goal= 0.0f;
+
+               verts->flags = 0;
+               VECCOPY ( verts->xold, verts->x );
+               VECCOPY ( verts->xconst, verts->x );
+               VECCOPY ( verts->txold, verts->x );
+               VecMulf ( verts->v, 0.0f );
+
+               verts->impulse_count = 0;
+               VECCOPY ( verts->impulse, tnull );
+       }
+       
+       // apply / set vertex groups
+       // has to be happen before springs are build!
+       cloth_apply_vgroup (clmd, dm);
+       
+       if ( !cloth_build_springs ( clmd, dm ) )
+       {
+               cloth_free_modifier ( ob, clmd );
+               modifier_setError ( & ( clmd->modifier ), "Can't build springs." );
+               printf("cloth_free_modifier cloth_build_springs\n");
+               return 0;
+       }
+       
+       // init our solver
+       if ( solvers [clmd->sim_parms->solver_type].init )
+               solvers [clmd->sim_parms->solver_type].init ( ob, clmd );
+
+       clmd->clothObject->tree = bvh_build_from_cloth ( clmd, clmd->coll_parms->epsilon );
+       
+       return 1;
+}
+
+
+static void cloth_from_mesh ( Object *ob, ClothModifierData *clmd, DerivedMesh *dm )
+{
+       unsigned int numverts = dm->getNumVerts ( dm );
+       unsigned int numfaces = dm->getNumFaces ( dm );
+       MFace *mface = CDDM_get_faces(dm);
+       unsigned int i = 0;
+
+       /* Allocate our vertices.
+       */
+       clmd->clothObject->numverts = numverts;
+       clmd->clothObject->verts = MEM_callocN ( sizeof ( ClothVertex ) * clmd->clothObject->numverts, "clothVertex" );
+       if ( clmd->clothObject->verts == NULL )
+       {
+               cloth_free_modifier ( ob, clmd );
+               modifier_setError ( & ( clmd->modifier ), "Out of memory on allocating clmd->clothObject->verts." );
+               printf("cloth_free_modifier clmd->clothObject->verts\n");
+               return;
+       }
+
+       // save face information
+       clmd->clothObject->numfaces = numfaces;
+       clmd->clothObject->mfaces = MEM_callocN ( sizeof ( MFace ) * clmd->clothObject->numfaces, "clothMFaces" );
+       if ( clmd->clothObject->mfaces == NULL )
+       {
+               cloth_free_modifier ( ob, clmd );
+               modifier_setError ( & ( clmd->modifier ), "Out of memory on allocating clmd->clothObject->mfaces." );
+               printf("cloth_free_modifier clmd->clothObject->mfaces\n");
+               return;
+       }
+       for ( i = 0; i < numfaces; i++ )
+               memcpy ( &clmd->clothObject->mfaces[i], &mface[i], sizeof ( MFace ) );
+
+       /* Free the springs since they can't be correct if the vertices
+       * changed.
+       */
+       if ( clmd->clothObject->springs != NULL )
+               MEM_freeN ( clmd->clothObject->springs );
+
+}
+
+/***************************************************************************************
+* SPRING NETWORK BUILDING IMPLEMENTATION BEGIN
+***************************************************************************************/
+
+// be carefull: implicit solver has to be resettet when using this one!
+// --> only for implicit handling of this spring!
+int cloth_add_spring ( ClothModifierData *clmd, unsigned int indexA, unsigned int indexB, float restlength, int spring_type)
+{
+       Cloth *cloth = clmd->clothObject;
+       ClothSpring *spring = NULL;
+       
+       if(cloth)
+       {
+               // TODO: look if this spring is already there
+               
+               spring = ( ClothSpring * ) MEM_callocN ( sizeof ( ClothSpring ), "cloth spring" );
+               
+               spring->ij = indexA;
+               spring->kl = indexB;
+               spring->restlen =  restlength;
+               spring->type = spring_type;
+               spring->flags = 0;
+               
+               cloth->numsprings++;
+       
+               BLI_linklist_append ( &cloth->springs, spring );
+               
+               return 1;
+       }
+       return 0;
+}
+
+int cloth_build_springs ( ClothModifierData *clmd, DerivedMesh *dm )
+{
+       Cloth *cloth = clmd->clothObject;
+       ClothSpring *spring = NULL, *tspring = NULL, *tspring2 = NULL;
+       unsigned int struct_springs = 0, shear_springs=0, bend_springs = 0;
+       unsigned int i = 0;
+       unsigned int numverts = dm->getNumVerts ( dm );
+       unsigned int numedges = dm->getNumEdges ( dm );
+       unsigned int numfaces = dm->getNumFaces ( dm );
+       MEdge *medge = CDDM_get_edges ( dm );
+       MFace *mface = CDDM_get_faces ( dm );
+       unsigned int index2 = 0; // our second vertex index
+       LinkNode **edgelist = NULL;
+       EdgeHash *edgehash = NULL;
+       LinkNode *search = NULL, *search2 = NULL;
+       float temp[3];
+       LinkNode *node = NULL, *node2 = NULL;
+       
+       // error handling
+       if ( numedges==0 )
+               return 0;
+
+       cloth->springs = NULL;
+
+       edgelist = MEM_callocN ( sizeof ( LinkNode * ) * numverts, "cloth_edgelist_alloc" );
+       for ( i = 0; i < numverts; i++ )
+       {
+               edgelist[i] = NULL;
+       }
+
+       if ( cloth->springs )
+               MEM_freeN ( cloth->springs );
+
+       // create spring network hash
+       edgehash = BLI_edgehash_new();
+
+       // structural springs
+       for ( i = 0; i < numedges; i++ )
+       {
+               spring = ( ClothSpring * ) MEM_callocN ( sizeof ( ClothSpring ), "cloth spring" );
+
+               if ( spring )
+               {
+                       spring->ij = medge[i].v1;
+                       spring->kl = medge[i].v2;
+                       VECSUB ( temp, cloth->verts[spring->kl].x, cloth->verts[spring->ij].x );
+                       spring->restlen =  sqrt ( INPR ( temp, temp ) );
+                       /*
+                       if(spring->restlen > 1.0)
+                       {
+                       printf("i: %d, L: %f\n", i, spring->restlen);
+                       printf("%d, x: %f, y: %f, z: %f\n", cloth->verts[spring->ij].x[0], cloth->verts[spring->ij].x[1], spring->ij, cloth->verts[spring->ij].x[2]);
+                       printf("%d, x: %f, y: %f, z: %f\n\n",spring->kl, cloth->verts[spring->kl].x[0], cloth->verts[spring->kl].x[1], cloth->verts[spring->kl].x[2]);
+               }
+                       */
+                       clmd->coll_parms->avg_spring_len += spring->restlen;
+                       spring->type = CLOTH_SPRING_TYPE_STRUCTURAL;
+                       spring->flags = 0;
+                       spring->stiffness = (cloth->verts[spring->kl].struct_stiff + cloth->verts[spring->ij].struct_stiff) / 2.0;
+                       struct_springs++;
+                       
+                       if(!i)
+                               node2 = BLI_linklist_append_fast ( &cloth->springs, spring );
+                       else
+                               node2 = BLI_linklist_append_fast ( &node->next, spring );
+                       node = node2;
+               }
+       }
+       
+       clmd->coll_parms->avg_spring_len /= struct_springs;
+       
+       // shear springs
+       for ( i = 0; i < numfaces; i++ )
+       {
+               spring = ( ClothSpring *) MEM_callocN ( sizeof ( ClothSpring ), "cloth spring" );
+
+               spring->ij = mface[i].v1;
+               spring->kl = mface[i].v3;
+               VECSUB ( temp, cloth->verts[spring->kl].x, cloth->verts[spring->ij].x );
+               spring->restlen =  sqrt ( INPR ( temp, temp ) );
+               spring->type = CLOTH_SPRING_TYPE_SHEAR;
+               spring->stiffness = (cloth->verts[spring->kl].shear_stiff + cloth->verts[spring->ij].shear_stiff) / 2.0;
+
+               BLI_linklist_append ( &edgelist[spring->ij], spring );
+               BLI_linklist_append ( &edgelist[spring->kl], spring );
+               shear_springs++;
+
+               node2 = BLI_linklist_append_fast ( &node->next, spring );
+               node = node2;
+
+               if ( mface[i].v4 )
+               {
+                       spring = ( ClothSpring * ) MEM_callocN ( sizeof ( ClothSpring ), "cloth spring" );
+
+                       spring->ij = mface[i].v2;
+                       spring->kl = mface[i].v4;
+                       VECSUB ( temp, cloth->verts[spring->kl].x, cloth->verts[spring->ij].x );
+                       spring->restlen =  sqrt ( INPR ( temp, temp ) );
+                       spring->type = CLOTH_SPRING_TYPE_SHEAR;
+                       spring->stiffness = (cloth->verts[spring->kl].shear_stiff + cloth->verts[spring->ij].shear_stiff) / 2.0;
+
+                       BLI_linklist_append ( &edgelist[spring->ij], spring );
+                       BLI_linklist_append ( &edgelist[spring->kl], spring );
+                       shear_springs++;
+
+                       node2 = BLI_linklist_append_fast ( &node->next, spring );
+                       node = node2;
+               }
+       }
+       
+       // bending springs
+       search2 = cloth->springs;
+       for ( i = struct_springs; i < struct_springs+shear_springs; i++ )
+       {
+               if ( !search2 )
+                       break;
+
+               tspring2 = search2->link;
+               search = edgelist[tspring2->kl];
+               while ( search )
+               {
+                       tspring = search->link;
+                       index2 = ( ( tspring->ij==tspring2->kl ) ? ( tspring->kl ) : ( tspring->ij ) );
+                       
+                       // check for existing spring
+                       // check also if startpoint is equal to endpoint
+                       if ( !BLI_edgehash_haskey ( edgehash, index2, tspring2->ij )
+                                                  && !BLI_edgehash_haskey ( edgehash, tspring2->ij, index2 )
+                                                  && ( index2!=tspring2->ij ) )
+                       {
+                               spring = ( ClothSpring * ) MEM_callocN ( sizeof ( ClothSpring ), "cloth spring" );
+
+                               spring->ij = tspring2->ij;
+                               spring->kl = index2;
+                               VECSUB ( temp, cloth->verts[index2].x, cloth->verts[tspring2->ij].x );
+                               spring->restlen =  sqrt ( INPR ( temp, temp ) );
+                               spring->type = CLOTH_SPRING_TYPE_BENDING;
+                               spring->stiffness = (cloth->verts[spring->kl].bend_stiff + cloth->verts[spring->ij].bend_stiff) / 2.0;
+                               BLI_edgehash_insert ( edgehash, spring->ij, index2, NULL );
+                               bend_springs++;
+
+                               node2 = BLI_linklist_append_fast ( &node->next, spring );
+                               node = node2;
+                       }
+                       search = search->next;
+               }
+               search2 = search2->next;
+       }
+       
+       cloth->numsprings = struct_springs + shear_springs + bend_springs;
+       
+       for ( i = 0; i < numverts; i++ )
+       {
+               BLI_linklist_free ( edgelist[i],NULL );
+       }
+       if ( edgelist )
+               MEM_freeN ( edgelist );
+       
+       BLI_edgehash_free ( edgehash, NULL );
+
+       return 1;
+
+} /* cloth_build_springs */
+/***************************************************************************************
+* SPRING NETWORK BUILDING IMPLEMENTATION END
+***************************************************************************************/
+
diff --git a/source/blender/blenkernel/intern/collision.c b/source/blender/blenkernel/intern/collision.c
new file mode 100644 (file)
index 0000000..5946b84
--- /dev/null
@@ -0,0 +1,1261 @@
+/*  collision.c      
+* 
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) Blender Foundation
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): none yet.
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "MEM_guardedalloc.h"
+/* types */
+#include "DNA_curve_types.h"
+#include "DNA_object_types.h"
+#include "DNA_object_force.h"
+#include "DNA_cloth_types.h"   
+#include "DNA_key_types.h"
+#include "DNA_mesh_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_edgehash.h"
+#include "BLI_linklist.h"
+#include "BKE_curve.h"
+#include "BKE_deform.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_displist.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+#include "BKE_mesh.h"
+#include "BKE_object.h"
+#include "BKE_cloth.h"
+#include "BKE_modifier.h"
+#include "BKE_utildefines.h"
+#include "BKE_DerivedMesh.h"
+#include "DNA_screen_types.h"
+#include "BSE_headerbuttons.h"
+#include "BIF_screen.h"
+#include "BIF_space.h"
+#include "mydevice.h"
+
+#include "Bullet-C-Api.h"
+
+/***********************************
+Collision modifier code start
+***********************************/
+
+/* step is limited from 0 (frame start position) to 1 (frame end position) */
+void collision_move_object(CollisionModifierData *collmd, float step, float prevstep)
+{
+       float tv[3] = {0,0,0};
+       unsigned int i = 0;
+       
+       for ( i = 0; i < collmd->numverts; i++ )
+       {
+               VECSUB(tv, collmd->xnew[i].co, collmd->x[i].co);
+               VECADDS(collmd->current_x[i].co, collmd->x[i].co, tv, prevstep);
+               VECADDS(collmd->current_xnew[i].co, collmd->x[i].co, tv, step);
+               VECSUB(collmd->current_v[i].co, collmd->current_xnew[i].co, collmd->current_x[i].co);
+       }
+}
+
+/* build bounding volume hierarchy from mverts (see kdop.c for whole BVH code) */
+BVH *bvh_build_from_mvert (MFace *mfaces, unsigned int numfaces, MVert *x, unsigned int numverts, float epsilon)
+{
+       BVH *bvh=NULL;
+       
+       bvh = MEM_callocN(sizeof(BVH), "BVH");
+       if (bvh == NULL) 
+       {
+               printf("bvh: Out of memory.\n");
+               return NULL;
+       }
+       
+       bvh->flags = 0;
+       bvh->leaf_tree = NULL;
+       bvh->leaf_root = NULL;
+       bvh->tree = NULL;
+
+       bvh->epsilon = epsilon;
+       bvh->numfaces = numfaces;
+       bvh->mfaces = mfaces;
+       
+       // we have no faces, we save seperate points
+       if(!mfaces)
+       {
+               bvh->numfaces = numverts;
+       }
+
+       bvh->numverts = numverts;
+       bvh->current_x = MEM_dupallocN(x);      
+       bvh->current_xold = MEM_dupallocN(x);   
+       
+       bvh_build(bvh);
+       
+       return bvh;
+}
+
+void bvh_update_from_mvert(BVH * bvh, MVert *x, unsigned int numverts, MVert *xnew, int moving)
+{
+       if(!bvh)
+               return;
+       
+       if(numverts!=bvh->numverts)
+               return;
+       
+       if(x)
+               memcpy(bvh->current_xold, x, sizeof(MVert) * numverts);
+       
+       if(xnew)
+               memcpy(bvh->current_x, xnew, sizeof(MVert) * numverts);
+       
+       bvh_update(bvh, moving);
+}
+
+/***********************************
+Collision modifier code end
+***********************************/
+
+/**
+ * gsl_poly_solve_cubic -
+ *
+ * copied from SOLVE_CUBIC.C --> GSL
+ */
+#define mySWAP(a,b) { float tmp = b ; b = a ; a = tmp ; }
+
+int gsl_poly_solve_cubic (float a, float b, float c, float *x0, float *x1, float *x2)
+{
+       float q = (a * a - 3 * b);
+       float r = (2 * a * a * a - 9 * a * b + 27 * c);
+
+       float Q = q / 9;
+       float R = r / 54;
+
+       float Q3 = Q * Q * Q;
+       float R2 = R * R;
+
+       float CR2 = 729 * r * r;
+       float CQ3 = 2916 * q * q * q;
+
+       if (R == 0 && Q == 0)
+       {
+               *x0 = - a / 3 ;
+               *x1 = - a / 3 ;
+               *x2 = - a / 3 ;
+               return 3 ;
+       }
+       else if (CR2 == CQ3) 
+       {
+         /* this test is actually R2 == Q3, written in a form suitable
+               for exact computation with integers */
+
+         /* Due to finite precision some float roots may be missed, and
+               considered to be a pair of complex roots z = x +/- epsilon i
+               close to the real axis. */
+
+               float sqrtQ = sqrtf (Q);
+
+               if (R > 0)
+               {
+                       *x0 = -2 * sqrtQ  - a / 3;
+                       *x1 = sqrtQ - a / 3;
+                       *x2 = sqrtQ - a / 3;
+               }
+               else
+               {
+                       *x0 = - sqrtQ  - a / 3;
+                       *x1 = - sqrtQ - a / 3;
+                       *x2 = 2 * sqrtQ - a / 3;
+               }
+               return 3 ;
+       }
+       else if (CR2 < CQ3) /* equivalent to R2 < Q3 */
+       {
+               float sqrtQ = sqrtf (Q);
+               float sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
+               float theta = acosf (R / sqrtQ3);
+               float norm = -2 * sqrtQ;
+               *x0 = norm * cosf (theta / 3) - a / 3;
+               *x1 = norm * cosf ((theta + 2.0 * M_PI) / 3) - a / 3;
+               *x2 = norm * cosf ((theta - 2.0 * M_PI) / 3) - a / 3;
+      
+               /* Sort *x0, *x1, *x2 into increasing order */
+
+               if (*x0 > *x1)
+                       mySWAP(*x0, *x1) ;
+      
+               if (*x1 > *x2)
+               {
+                       mySWAP(*x1, *x2) ;
+          
+                       if (*x0 > *x1)
+                               mySWAP(*x0, *x1) ;
+               }
+      
+               return 3;
+       }
+       else
+       {
+               float sgnR = (R >= 0 ? 1 : -1);
+               float A = -sgnR * powf (fabs (R) + sqrtf (R2 - Q3), 1.0/3.0);
+               float B = Q / A ;
+               *x0 = A + B - a / 3;
+               return 1;
+       }
+}
+
+
+/**
+ * gsl_poly_solve_quadratic
+ *
+ * copied from GSL
+ */
+int gsl_poly_solve_quadratic (float a, float b, float c,  float *x0, float *x1)
+{
+       float disc = b * b - 4 * a * c;
+
+       if (disc > 0)
+       {
+               if (b == 0)
+               {
+                       float r = fabs (0.5 * sqrtf (disc) / a);
+                       *x0 = -r;
+                       *x1 =  r;
+               }
+               else
+               {
+                       float sgnb = (b > 0 ? 1 : -1);
+                       float temp = -0.5 * (b + sgnb * sqrtf (disc));
+                       float r1 = temp / a ;
+                       float r2 = c / temp ;
+
+                       if (r1 < r2) 
+                       {
+                               *x0 = r1 ;
+                               *x1 = r2 ;
+                       } 
+                       else 
+                       {
+                               *x0 = r2 ;
+                               *x1 = r1 ;
+                       }
+               }
+               return 2;
+       }
+       else if (disc == 0) 
+       {
+               *x0 = -0.5 * b / a ;
+               *x1 = -0.5 * b / a ;
+               return 2 ;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+
+
+/*
+ * See Bridson et al. "Robust Treatment of Collision, Contact and Friction for Cloth Animation"
+ *     page 4, left column
+ */
+
+int cloth_get_collision_time(float a[3], float b[3], float c[3], float d[3], float e[3], float f[3], float solution[3]) 
+{
+       int num_sols = 0;
+       
+       float g = -a[2] * c[1] * e[0] + a[1] * c[2] * e[0] +
+                       a[2] * c[0] * e[1] - a[0] * c[2] * e[1] -
+                       a[1] * c[0] * e[2] + a[0] * c[1] * e[2];
+
+       float h = -b[2] * c[1] * e[0] + b[1] * c[2] * e[0] - a[2] * d[1] * e[0] +
+                       a[1] * d[2] * e[0] + b[2] * c[0] * e[1] - b[0] * c[2] * e[1] +
+                       a[2] * d[0] * e[1] - a[0] * d[2] * e[1] - b[1] * c[0] * e[2] +
+                       b[0] * c[1] * e[2] - a[1] * d[0] * e[2] + a[0] * d[1] * e[2] -
+                       a[2] * c[1] * f[0] + a[1] * c[2] * f[0] + a[2] * c[0] * f[1] -
+                       a[0] * c[2] * f[1] - a[1] * c[0] * f[2] + a[0] * c[1] * f[2];
+
+       float i = -b[2] * d[1] * e[0] + b[1] * d[2] * e[0] +
+                       b[2] * d[0] * e[1] - b[0] * d[2] * e[1] -
+                       b[1] * d[0] * e[2] + b[0] * d[1] * e[2] -
+                       b[2] * c[1] * f[0] + b[1] * c[2] * f[0] -
+                       a[2] * d[1] * f[0] + a[1] * d[2] * f[0] +
+                       b[2] * c[0] * f[1] - b[0] * c[2] * f[1] + 
+                       a[2] * d[0] * f[1] - a[0] * d[2] * f[1] -
+                       b[1] * c[0] * f[2] + b[0] * c[1] * f[2] -
+                       a[1] * d[0] * f[2] + a[0] * d[1] * f[2];
+
+       float j = -b[2] * d[1] * f[0] + b[1] * d[2] * f[0] +
+                       b[2] * d[0] * f[1] - b[0] * d[2] * f[1] -
+                       b[1] * d[0] * f[2] + b[0] * d[1] * f[2];
+
+       // Solve cubic equation to determine times t1, t2, t3, when the collision will occur.
+       if(ABS(j) > ALMOST_ZERO)
+       {
+               i /= j;
+               h /= j;
+               g /= j;
+               
+               num_sols = gsl_poly_solve_cubic(i, h, g, &solution[0], &solution[1], &solution[2]);
+       }
+       else if(ABS(i) > ALMOST_ZERO)
+       {       
+               num_sols = gsl_poly_solve_quadratic(i, h, g, &solution[0], &solution[1]);
+               solution[2] = -1.0;
+       }
+       else if(ABS(h) > ALMOST_ZERO)
+       {
+               solution[0] = -g / h;
+               solution[1] = solution[2] = -1.0;
+               num_sols = 1;
+       }
+       else if(ABS(g) > ALMOST_ZERO)
+       {
+               solution[0] = 0;
+               solution[1] = solution[2] = -1.0;
+               num_sols = 1;
+       }
+
+       // Discard negative solutions
+       if ((num_sols >= 1) && (solution[0] < 0)) 
+       {
+               --num_sols;
+               solution[0] = solution[num_sols];
+       }
+       if ((num_sols >= 2) && (solution[1] < 0)) 
+       {
+               --num_sols;
+               solution[1] = solution[num_sols];
+       }
+       if ((num_sols == 3) && (solution[2] < 0)) 
+       {
+               --num_sols;
+       }
+
+       // Sort
+       if (num_sols == 2) 
+       {
+               if (solution[0] > solution[1]) 
+               {
+                       double tmp = solution[0];
+                       solution[0] = solution[1];
+                       solution[1] = tmp;
+               }
+       }
+       else if (num_sols == 3) 
+       {
+
+               // Bubblesort
+               if (solution[0] > solution[1]) {
+                       double tmp = solution[0]; solution[0] = solution[1]; solution[1] = tmp;
+               }
+               if (solution[1] > solution[2]) {
+                       double tmp = solution[1]; solution[1] = solution[2]; solution[2] = tmp;
+               }
+               if (solution[0] > solution[1]) {
+                       double tmp = solution[0]; solution[0] = solution[1]; solution[1] = tmp;
+               }
+       }
+
+       return num_sols;
+}
+
+// w3 is not perfect
+void cloth_compute_barycentric (float pv[3], float p1[3], float p2[3], float p3[3], float *w1, float *w2, float *w3)
+{
+       double  tempV1[3], tempV2[3], tempV4[3];
+       double  a,b,c,d,e,f;
+
+       VECSUB (tempV1, p1, p3);        
+       VECSUB (tempV2, p2, p3);        
+       VECSUB (tempV4, pv, p3);        
+       
+       a = INPR (tempV1, tempV1);      
+       b = INPR (tempV1, tempV2);      
+       c = INPR (tempV2, tempV2);      
+       e = INPR (tempV1, tempV4);      
+       f = INPR (tempV2, tempV4);      
+       
+       d = (a * c - b * b);
+       
+       if (ABS(d) < ALMOST_ZERO) {
+               *w1 = *w2 = *w3 = 1.0 / 3.0;
+               return;
+       }
+       
+       w1[0] = (float)((e * c - b * f) / d);
+       
+       if(w1[0] < 0)
+               w1[0] = 0;
+       
+       w2[0] = (float)((f - b * (double)w1[0]) / c);
+       
+       if(w2[0] < 0)
+               w2[0] = 0;
+       
+       w3[0] = 1.0f - w1[0] - w2[0];
+}
+
+DO_INLINE void interpolateOnTriangle(float to[3], float v1[3], float v2[3], float v3[3], double w1, double w2, double w3) 
+{
+       to[0] = to[1] = to[2] = 0;
+       VECADDMUL(to, v1, w1);
+       VECADDMUL(to, v2, w2);
+       VECADDMUL(to, v3, w3);
+}
+
+// unused in the moment, has some bug in
+DO_INLINE void calculateFrictionImpulse(float to[3], float vrel[3], float normal[3], double normalVelocity,
+                                       double frictionConstant, double delta_V_n) 
+{
+       float vrel_t_pre[3];
+       float vrel_t[3];
+       VECSUBS(vrel_t_pre, vrel, normal, normalVelocity);
+       VECCOPY(to, vrel_t_pre);
+       VecMulf(to, MAX2(1.0f - frictionConstant * delta_V_n / INPR(vrel_t_pre,vrel_t_pre), 0.0f));
+}
+
+int cloth_collision_response_static(ClothModifierData *clmd, CollisionModifierData *collmd)
+{
+       unsigned int i = 0;
+       int result = 0;
+       LinkNode *search = NULL;
+       CollPair *collpair = NULL;
+       Cloth *cloth1;
+       float w1, w2, w3, u1, u2, u3;
+       float v1[3], v2[3], relativeVelocity[3];
+       float magrelVel;
+       
+       cloth1 = clmd->clothObject;
+
+       search = clmd->coll_parms->collision_list;
+       
+       while(search)
+       {
+               collpair = search->link;
+               
+               // compute barycentric coordinates for both collision points
+               cloth_compute_barycentric(collpair->pa,
+                                         cloth1->verts[collpair->ap1].txold,
+       cloth1->verts[collpair->ap2].txold,
+       cloth1->verts[collpair->ap3].txold, 
+       &w1, &w2, &w3);
+               
+               // was: txold
+               cloth_compute_barycentric(collpair->pb,
+                                         collmd->current_x[collpair->bp1].co,
+       collmd->current_x[collpair->bp2].co,
+       collmd->current_x[collpair->bp3].co,
+       &u1, &u2, &u3);
+       
+               // Calculate relative "velocity".
+               interpolateOnTriangle(v1, cloth1->verts[collpair->ap1].tv, cloth1->verts[collpair->ap2].tv, cloth1->verts[collpair->ap3].tv, w1, w2, w3);
+               
+               interpolateOnTriangle(v2, collmd->current_v[collpair->bp1].co, collmd->current_v[collpair->bp2].co, collmd->current_v[collpair->bp3].co, u1, u2, u3);
+               
+               VECSUB(relativeVelocity, v1, v2);
+                       
+               // Calculate the normal component of the relative velocity (actually only the magnitude - the direction is stored in 'normal').
+               magrelVel = INPR(relativeVelocity, collpair->normal);
+               
+               // printf("magrelVel: %f\n", magrelVel);
+                               
+               // Calculate masses of points.
+               
+               // If v_n_mag < 0 the edges are approaching each other.
+               if(magrelVel < -ALMOST_ZERO) 
+               {
+                       // Calculate Impulse magnitude to stop all motion in normal direction.
+                       // const double I_mag = v_n_mag / (1/m1 + 1/m2);
+                       float magnitude_i = magrelVel / 2.0f; // TODO implement masses
+                       float tangential[3], magtangent, magnormal, collvel[3];
+                       float vrel_t_pre[3];
+                       float vrel_t[3];
+                       double impulse;
+                       float epsilon = clmd->coll_parms->epsilon;
+                       float overlap = (epsilon + ALMOST_ZERO-collpair->distance);
+                       
+                       // calculateFrictionImpulse(tangential, relativeVelocity, collpair->normal, magrelVel, clmd->coll_parms->friction*0.01, magrelVel);
+                       
+                       // magtangent = INPR(tangential, tangential);
+                       
+                       // Apply friction impulse.
+                       if (magtangent < -ALMOST_ZERO) 
+                       {
+                               
+                               // printf("friction applied: %f\n", magtangent);
+                               // TODO check original code 
+                               /*
+                               VECSUB(cloth1->verts[face1->v1].tv, cloth1->verts[face1->v1].tv,tangential);
+                               VECSUB(cloth1->verts[face1->v1].tv, cloth1->verts[face1->v2].tv,tangential);
+                               VECSUB(cloth1->verts[face1->v1].tv, cloth1->verts[face1->v3].tv,tangential);
+                               VECSUB(cloth1->verts[face1->v1].tv, cloth1->verts[face1->v4].tv,tangential);
+                               */
+                       }
+                       
+
+                       impulse = -2.0f * magrelVel / ( 1.0 + w1*w1 + w2*w2 + w3*w3);
+                       
+                       // printf("impulse: %f\n", impulse);
+                       
+                       VECADDMUL(cloth1->verts[collpair->ap1].impulse, collpair->normal, w1 * impulse); 
+                       cloth1->verts[collpair->ap1].impulse_count++;
+                       
+                       VECADDMUL(cloth1->verts[collpair->ap2].impulse, collpair->normal, w2 * impulse); 
+                       cloth1->verts[collpair->ap2].impulse_count++;
+                       
+                       VECADDMUL(cloth1->verts[collpair->ap3].impulse, collpair->normal, w3 * impulse); 
+                       cloth1->verts[collpair->ap3].impulse_count++;
+                       
+                       result = 1;
+                       
+                       /*
+                       if (overlap > ALMOST_ZERO) {
+                       double I_mag  = overlap * 0.1;
+                               
+                       impulse = -I_mag / ( 1.0 + w1*w1 + w2*w2 + w3*w3);
+                               
+                       VECADDMUL(cloth1->verts[collpair->ap1].impulse, collpair->normal, w1 * impulse); 
+                       cloth1->verts[collpair->ap1].impulse_count++;
+                                                       
+                       VECADDMUL(cloth1->verts[collpair->ap2].impulse, collpair->normal, w2 * impulse); 
+                       cloth1->verts[collpair->ap2].impulse_count++;
+                       
+                       VECADDMUL(cloth1->verts[collpair->ap3].impulse, collpair->normal, w3 * impulse); 
+                       cloth1->verts[collpair->ap3].impulse_count++;
+               }
+                       */
+               
+                       // printf("magnitude_i: %f\n", magnitude_i); // negative before collision in my case
+                       
+                       // Apply the impulse and increase impulse counters.
+
+                       /*              
+                       // calculateFrictionImpulse(tangential, collvel, collpair->normal, magtangent, clmd->coll_parms->friction*0.01, magtangent);
+                       VECSUBS(vrel_t_pre, collvel, collpair->normal, magnormal);
+                       // VecMulf(vrel_t_pre, clmd->coll_parms->friction*0.01f/INPR(vrel_t_pre,vrel_t_pre));
+                       magtangent = Normalize(vrel_t_pre);
+                       VecMulf(vrel_t_pre, MIN2(clmd->coll_parms->friction*0.01f*magnormal,magtangent));
+                       
+                       VECSUB(cloth1->verts[face1->v1].tv, cloth1->verts[face1->v1].tv,vrel_t_pre);
+                       */
+                       
+                       
+                       
+               }
+               
+               search = search->next;
+       }
+       
+               
+       return result;
+}
+
+int cloth_collision_response_moving_tris(ClothModifierData *clmd, ClothModifierData *coll_clmd)
+{
+       return 1;
+}
+
+
+int cloth_collision_response_moving_edges(ClothModifierData *clmd, ClothModifierData *coll_clmd)
+{
+       return 1;
+}
+
+void cloth_collision_static(ClothModifierData *clmd, CollisionModifierData *collmd, CollisionTree *tree1, CollisionTree *tree2)
+{
+       CollPair *collpair = NULL;
+       Cloth *cloth1=NULL;
+       MFace *face1=NULL, *face2=NULL;
+       ClothVertex *verts1=NULL;
+       double distance = 0;
+       float epsilon = clmd->coll_parms->epsilon;
+       unsigned int i = 0;
+
+       for(i = 0; i < 4; i++)
+       {
+               collpair = (CollPair *)MEM_callocN(sizeof(CollPair), "cloth coll pair");                
+               
+               cloth1 = clmd->clothObject;
+               
+               verts1 = cloth1->verts;
+       
+               face1 = &(cloth1->mfaces[tree1->tri_index]);
+               face2 = &(collmd->mfaces[tree2->tri_index]);
+               
+               // check all possible pairs of triangles
+               if(i == 0)
+               {
+                       collpair->ap1 = face1->v1;
+                       collpair->ap2 = face1->v2;
+                       collpair->ap3 = face1->v3;
+                       
+                       collpair->bp1 = face2->v1;
+                       collpair->bp2 = face2->v2;
+                       collpair->bp3 = face2->v3;
+                       
+               }
+               
+               if(i == 1)
+               {
+                       if(face1->v4)
+                       {
+                               collpair->ap1 = face1->v3;
+                               collpair->ap2 = face1->v4;
+                               collpair->ap3 = face1->v1;
+                               
+                               collpair->bp1 = face2->v1;
+                               collpair->bp2 = face2->v2;
+                               collpair->bp3 = face2->v3;
+                       }
+                       else
+                               i++;
+               }
+               
+               if(i == 2)
+               {
+                       if(face2->v4)
+                       {
+                               collpair->ap1 = face1->v1;
+                               collpair->ap2 = face1->v2;
+                               collpair->ap3 = face1->v3;
+                               
+                               collpair->bp1 = face2->v3;
+                               collpair->bp2 = face2->v4;
+                               collpair->bp3 = face2->v1;
+                       }
+                       else
+                               i+=2;
+               }
+               
+               if(i == 3)
+               {
+                       if((face1->v4)&&(face2->v4))
+                       {
+                               collpair->ap1 = face1->v3;
+                               collpair->ap2 = face1->v4;
+                               collpair->ap3 = face1->v1;
+                               
+                               collpair->bp1 = face2->v3;
+                               collpair->bp2 = face2->v4;
+                               collpair->bp3 = face2->v1;
+                       }
+                       else
+                               i++;
+               }
+               
+               // calc SIPcode (?)
+               
+               if(i < 4)
+               {
+                       // calc distance + normal       
+                       distance = plNearestPoints(
+                                       verts1[collpair->ap1].txold, verts1[collpair->ap2].txold, verts1[collpair->ap3].txold, collmd->current_x[collpair->bp1].co, collmd->current_x[collpair->bp2].co, collmd->current_x[collpair->bp3].co, collpair->pa,collpair->pb,collpair->vector);
+                       
+                       if (distance <= (epsilon + ALMOST_ZERO))
+                       {
+                               // printf("dist: %f\n", (float)distance);
+                               
+                               // collpair->face1 = tree1->tri_index;
+                               // collpair->face2 = tree2->tri_index;
+                               
+                               VECCOPY(collpair->normal, collpair->vector);
+                               Normalize(collpair->normal);
+                               
+                               collpair->distance = distance;
+                               BLI_linklist_prepend(&clmd->coll_parms->collision_list, collpair);
+                               
+                       }
+                       else
+                       {
+                               MEM_freeN(collpair);
+                       }
+               }
+               else
+               {
+                       MEM_freeN(collpair);
+               }
+       }
+}
+
+int cloth_are_edges_adjacent(ClothModifierData *clmd, ClothModifierData *coll_clmd, EdgeCollPair *edgecollpair)
+{
+       Cloth *cloth1 = NULL, *cloth2 = NULL;
+       ClothVertex *verts1 = NULL, *verts2 = NULL;
+       float temp[3];
+        
+       cloth1 = clmd->clothObject;
+       cloth2 = coll_clmd->clothObject;
+       
+       verts1 = cloth1->verts;
+       verts2 = cloth2->verts;
+       
+       VECSUB(temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p21].xold);
+       if(ABS(INPR(temp, temp)) < ALMOST_ZERO)
+               return 1;
+       
+       VECSUB(temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p22].xold);
+       if(ABS(INPR(temp, temp)) < ALMOST_ZERO)
+               return 1;
+       
+       VECSUB(temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p21].xold);
+       if(ABS(INPR(temp, temp)) < ALMOST_ZERO)
+               return 1;
+       
+       VECSUB(temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p22].xold);
+       if(ABS(INPR(temp, temp)) < ALMOST_ZERO)
+               return 1;
+               
+       return 0;
+}
+
+void cloth_collision_moving_edges(ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2)
+{
+       EdgeCollPair edgecollpair;
+       Cloth *cloth1=NULL, *cloth2=NULL;
+       MFace *face1=NULL, *face2=NULL;
+       ClothVertex *verts1=NULL, *verts2=NULL;
+       double distance = 0;
+       float epsilon = clmd->coll_parms->epsilon;
+       unsigned int i = 0, j = 0, k = 0;
+       int numsolutions = 0;
+       float a[3], b[3], c[3], d[3], e[3], f[3], solution[3];
+       
+       cloth1 = clmd->clothObject;
+       cloth2 = coll_clmd->clothObject;
+       
+       verts1 = cloth1->verts;
+       verts2 = cloth2->verts;
+
+       face1 = &(cloth1->mfaces[tree1->tri_index]);
+       face2 = &(cloth2->mfaces[tree2->tri_index]);
+       
+       for( i = 0; i < 5; i++)
+       {
+               if(i == 0) 
+               {
+                       edgecollpair.p11 = face1->v1;
+                       edgecollpair.p12 = face1->v2;
+               }
+               else if(i == 1) 
+               {
+                       edgecollpair.p11 = face1->v2;
+                       edgecollpair.p12 = face1->v3;
+               }
+               else if(i == 2) 
+               {
+                       if(face1->v4) 
+                       {
+                               edgecollpair.p11 = face1->v3;
+                               edgecollpair.p12 = face1->v4;
+                       }
+                       else 
+                       {
+                               edgecollpair.p11 = face1->v3;
+                               edgecollpair.p12 = face1->v1;
+                               i+=5; // get out of here after this edge pair is handled
+                       }
+               }
+               else if(i == 3) 
+               {
+                       if(face1->v4) 
+                       {
+                               edgecollpair.p11 = face1->v4;
+                               edgecollpair.p12 = face1->v1;
+                       }       
+                       else
+                               continue;
+               }
+               else
+               {
+                       edgecollpair.p11 = face1->v3;
+                       edgecollpair.p12 = face1->v1;
+               }
+
+               
+               for( j = 0; j < 5; j++)
+               {
+                       if(j == 0)
+                       {
+                               edgecollpair.p21 = face2->v1;
+                               edgecollpair.p22 = face2->v2;
+                       }
+                       else if(j == 1)
+                       {
+                               edgecollpair.p21 = face2->v2;
+                               edgecollpair.p22 = face2->v3;
+                       }
+                       else if(j == 2)
+                       {
+                               if(face2->v4) 
+                               {
+                                       edgecollpair.p21 = face2->v3;
+                                       edgecollpair.p22 = face2->v4;
+                               }
+                               else 
+                               {
+                                       edgecollpair.p21 = face2->v3;
+                                       edgecollpair.p22 = face2->v1;
+                               }
+                       }
+                       else if(j == 3)
+                       {
+                               if(face2->v4) 
+                               {
+                                       edgecollpair.p21 = face2->v4;
+                                       edgecollpair.p22 = face2->v1;
+                               }
+                               else
+                                       continue;
+                       }
+                       else
+                       {
+                               edgecollpair.p21 = face2->v3;
+                               edgecollpair.p22 = face2->v1;
+                       }
+                       
+                       
+                       if(!cloth_are_edges_adjacent(clmd, coll_clmd, &edgecollpair))
+                       {
+                               VECSUB(a, verts1[edgecollpair.p12].xold, verts1[edgecollpair.p11].xold);
+                               VECSUB(b, verts1[edgecollpair.p12].v, verts1[edgecollpair.p11].v);
+                               VECSUB(c, verts1[edgecollpair.p21].xold, verts1[edgecollpair.p11].xold);
+                               VECSUB(d, verts1[edgecollpair.p21].v, verts1[edgecollpair.p11].v);
+                               VECSUB(e, verts2[edgecollpair.p22].xold, verts1[edgecollpair.p11].xold);
+                               VECSUB(f, verts2[edgecollpair.p22].v, verts1[edgecollpair.p11].v);
+                               
+                               numsolutions = cloth_get_collision_time(a, b, c, d, e, f, solution);
+                               
+                               for (k = 0; k < numsolutions; k++) 
+                               {                                                               
+                                       if ((solution[k] >= 0.0) && (solution[k] <= 1.0)) 
+                                       {
+                                               float out_collisionTime = solution[k];
+                                               
+                                               // TODO: check for collisions 
+                                               
+                                               // TODO: put into (edge) collision list
+                                               
+                                               // printf("Moving edge found!\n");
+                                       }
+                               }
+                       }
+               }
+       }               
+}
+
+void cloth_collision_moving_tris(ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2)
+{
+       CollPair collpair;
+       Cloth *cloth1=NULL, *cloth2=NULL;
+       MFace *face1=NULL, *face2=NULL;
+       ClothVertex *verts1=NULL, *verts2=NULL;
+       double distance = 0;
+       float epsilon = clmd->coll_parms->epsilon;
+       unsigned int i = 0, j = 0, k = 0;
+       int numsolutions = 0;
+       float a[3], b[3], c[3], d[3], e[3], f[3], solution[3];
+
+       for(i = 0; i < 2; i++)
+       {               
+               cloth1 = clmd->clothObject;
+               cloth2 = coll_clmd->clothObject;
+               
+               verts1 = cloth1->verts;
+               verts2 = cloth2->verts;
+       
+               face1 = &(cloth1->mfaces[tree1->tri_index]);
+               face2 = &(cloth2->mfaces[tree2->tri_index]);
+               
+               // check all possible pairs of triangles
+               if(i == 0)
+               {
+                       collpair.ap1 = face1->v1;
+                       collpair.ap2 = face1->v2;
+                       collpair.ap3 = face1->v3;
+                       
+                       collpair.pointsb[0] = face2->v1;
+                       collpair.pointsb[1] = face2->v2;
+                       collpair.pointsb[2] = face2->v3;
+                       collpair.pointsb[3] = face2->v4;
+               }
+               
+               if(i == 1)
+               {
+                       if(face1->v4)
+                       {
+                               collpair.ap1 = face1->v3;
+                               collpair.ap2 = face1->v4;
+                               collpair.ap3 = face1->v1;
+                               
+                               collpair.pointsb[0] = face2->v1;
+                               collpair.pointsb[1] = face2->v2;
+                               collpair.pointsb[2] = face2->v3;
+                               collpair.pointsb[3] = face2->v4;
+                       }
+                       else
+                               i++;
+               }
+               
+               // calc SIPcode (?)
+               
+               if(i < 2)
+               {
+                       VECSUB(a, verts1[collpair.ap2].xold, verts1[collpair.ap1].xold);
+                       VECSUB(b, verts1[collpair.ap2].v, verts1[collpair.ap1].v);
+                       VECSUB(c, verts1[collpair.ap3].xold, verts1[collpair.ap1].xold);
+                       VECSUB(d, verts1[collpair.ap3].v, verts1[collpair.ap1].v);
+                               
+                       for(j = 0; j < 4; j++)
+                       {                                       
+                               if((j==3) && !(face2->v4))
+                                       break;
+                               
+                               VECSUB(e, verts2[collpair.pointsb[j]].xold, verts1[collpair.ap1].xold);
+                               VECSUB(f, verts2[collpair.pointsb[j]].v, verts1[collpair.ap1].v);
+                               
+                               numsolutions = cloth_get_collision_time(a, b, c, d, e, f, solution);
+                               
+                               for (k = 0; k < numsolutions; k++) 
+                               {                                                               
+                                       if ((solution[k] >= 0.0) && (solution[k] <= 1.0)) 
+                                       {
+                                               float out_collisionTime = solution[k];
+                                               
+                                               // TODO: check for collisions 
+                                               
+                                               // TODO: put into (point-face) collision list
+                                               
+                                               // printf("Moving found!\n");
+                                               
+                                       }
+                               }
+                               
+                               // TODO: check borders for collisions
+                       }
+                       
+               }
+       }
+}
+
+void cloth_collision_moving(ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2)
+{
+       // TODO: check for adjacent
+       cloth_collision_moving_edges(clmd, coll_clmd, tree1, tree2);
+       
+       cloth_collision_moving_tris(clmd, coll_clmd, tree1, tree2);
+       cloth_collision_moving_tris(coll_clmd, clmd, tree2, tree1);
+}
+
+// cloth - object collisions
+int cloth_bvh_objcollision(ClothModifierData * clmd, float step, float dt)
+{
+       Base *base=NULL;
+       CollisionModifierData *collmd=NULL;
+       Cloth *cloth=NULL;
+       Object *coll_ob=NULL;
+       BVH *cloth_bvh=NULL;
+       unsigned int i=0, j = 0, numfaces = 0, numverts = 0;
+       unsigned int result = 0, ic = 0, rounds = 0; // result counts applied collisions; ic is for debug output; 
+       ClothVertex *verts = NULL;
+       float tnull[3] = {0,0,0};
+       int ret = 0;
+       ClothModifierData *tclmd;
+
+       if ((clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ) || !(((Cloth *)clmd->clothObject)->tree))
+       {
+               return 0;
+       }
+       
+       cloth = clmd->clothObject;
+       verts = cloth->verts;
+       cloth_bvh = (BVH *) cloth->tree;
+       numfaces = clmd->clothObject->numfaces;
+       numverts = clmd->clothObject->numverts;
+       
+       ////////////////////////////////////////////////////////////
+       // static collisions
+       ////////////////////////////////////////////////////////////
+
+       // update cloth bvh
+       bvh_update_from_cloth(clmd, 0); // 0 means STATIC, 1 means MOVING (see later in this function)
+       
+       do
+       {
+               result = 0;
+               ic = 0;
+               clmd->coll_parms->collision_list = NULL; 
+               
+               // check all collision objects
+               for (base = G.scene->base.first; base; base = base->next)
+               {
+                       coll_ob = base->object;
+                       collmd = (CollisionModifierData *) modifiers_findByType (coll_ob, eModifierType_Collision);
+                       
+                       if (!collmd)
+                               continue;
+                       
+                       tclmd = (ClothModifierData *) modifiers_findByType (coll_ob, eModifierType_Cloth);
+                       if(tclmd == clmd)
+                               continue;
+                       
+                       if (collmd->tree) 
+                       {
+                               BVH *coll_bvh = collmd->tree;
+                               
+                               collision_move_object(collmd, step + dt, step);
+                                       
+                               bvh_traverse(clmd, collmd, cloth_bvh->root, coll_bvh->root, step, cloth_collision_static);
+                       }
+                       else
+                               printf ("cloth_bvh_objcollision: found a collision object with clothObject or collData NULL.\n");
+               
+               
+                       // process all collisions (calculate impulses, TODO: also repulses if distance too short)
+                       result = 1;
+                       for(j = 0; j < 10; j++) // 10 is just a value that ensures convergence
+                       {
+                               result = 0;
+                               
+                               if (collmd->tree) 
+                                       result += cloth_collision_response_static(clmd, collmd);
+                       
+                               
+                               // apply impulses in parallel
+                               ic=0;
+                               for(i = 0; i < numverts; i++)
+                               {
+                                       // calculate "velocities" (just xnew = xold + v; no dt in v)
+                                       if(verts[i].impulse_count)
+                                       {
+                                               VECADDMUL(verts[i].tv, verts[i].impulse, 1.0f / verts[i].impulse_count);
+                                               VECCOPY(verts[i].impulse, tnull);
+                                               verts[i].impulse_count = 0;
+                                       
+                                               ic++;
+                                               ret++;
+                                       }
+                               }
+                       }
+                       
+                       // free collision list
+                       if(clmd->coll_parms->collision_list)
+                       {
+                               LinkNode *search = clmd->coll_parms->collision_list;
+                               while(search)
+                               {
+                                       CollPair *coll_pair = search->link;
+                                                       
+                                       MEM_freeN(coll_pair);
+                                       search = search->next;
+                               }
+                               BLI_linklist_free(clmd->coll_parms->collision_list,NULL);
+                       
+                               clmd->coll_parms->collision_list = NULL;
+                       }
+               }
+               rounds++;
+       }
+       while(result && (clmd->coll_parms->loop_count>rounds));
+       
+       ////////////////////////////////////////////////////////////
+       // update positions
+       // this is needed for bvh_calc_DOP_hull_moving() [kdop.c]
+       ////////////////////////////////////////////////////////////
+       
+       // verts come from clmd
+       for(i = 0; i < numverts; i++)
+       {
+               if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) 
+               {                       
+                       if(verts [i].goal >= SOFTGOALSNAP)
+                       {
+                               continue;
+                       }
+               }
+               
+               VECADD(verts[i].tx, verts[i].txold, verts[i].tv);
+       }
+       ////////////////////////////////////////////////////////////
+
+       ////////////////////////////////////////////////////////////
+       // moving collisions
+       //
+       // response code is just missing itm 
+       ////////////////////////////////////////////////////////////
+
+       /*
+       // update cloth bvh
+       bvh_update_from_cloth(clmd, 1);  // 0 means STATIC, 1 means MOVING 
+       
+       // update moving bvh for collision object once
+       for (base = G.scene->base.first; base; base = base->next)
+       {
+               
+       coll_ob = base->object;
+       coll_clmd = (ClothModifierData *) modifiers_findByType (coll_ob, eModifierType_Cloth);
+       if (!coll_clmd)
+       continue;
+               
+       if(!coll_clmd->clothObject)
+       continue;
+               
+                               // if collision object go on
+       if (coll_clmd->clothObject && coll_clmd->clothObject->tree) 
+       {
+       BVH *coll_bvh = coll_clmd->clothObject->tree;
+                       
+       bvh_update_from_cloth(coll_clmd, 1);  // 0 means STATIC, 1 means MOVING         
+}
+}
+       
+       
+       do
+       {
+       result = 0;
+       ic = 0;
+       clmd->coll_parms->collision_list = NULL; 
+               
+               // check all collision objects
+       for (base = G.scene->base.first; base; base = base->next)
+       {
+       coll_ob = base->object;
+       coll_clmd = (ClothModifierData *) modifiers_findByType (coll_ob, eModifierType_Cloth);
+                       
+       if (!coll_clmd)
+       continue;
+                       
+                       // if collision object go on
+       if (coll_clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ)
+       {
+       if (coll_clmd->clothObject && coll_clmd->clothObject->tree) 
+       {
+       BVH *coll_bvh = coll_clmd->clothObject->tree;
+                                       
+       bvh_traverse(clmd, coll_clmd, cloth_bvh->root, coll_bvh->root, step, cloth_collision_moving);
+}
+       else
+       printf ("cloth_bvh_objcollision: found a collision object with clothObject or collData NULL.\n");
+}
+}
+               
+               // process all collisions (calculate impulses, TODO: also repulses if distance too short)
+       result = 1;
+       for(j = 0; j < 10; j++) // 10 is just a value that ensures convergence
+       {
+       result = 0;
+                               
+                       // handle all collision objects
+       for (base = G.scene->base.first; base; base = base->next)
+       {
+                       
+       coll_ob = base->object;
+       coll_clmd = (ClothModifierData *) modifiers_findByType (coll_ob, eModifierType_Cloth);
+                                               
+       if (!coll_clmd)
+       continue;
+                               
+                               // if collision object go on
+       if (coll_clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ)
+       {
+       if (coll_clmd->clothObject) 
+       result += cloth_collision_response_moving_tris(clmd, coll_clmd);
+       else
+       printf ("cloth_bvh_objcollision: found a collision object with clothObject or collData NULL.\n");
+}
+}
+                                               
+                       // apply impulses in parallel
+       ic=0;
+       for(i = 0; i < numverts; i++)
+       {
+                               // calculate "velocities" (just xnew = xold + v; no dt in v)
+       if(verts[i].impulse_count)
+       {
+       VECADDMUL(verts[i].tv, verts[i].impulse, 1.0f / verts[i].impulse_count);
+       VECCOPY(verts[i].impulse, tnull);
+       verts[i].impulse_count = 0;
+                                                       
+       ic++;
+       ret++;
+}
+}
+}
+               
+               
+               // verts come from clmd
+       for(i = 0; i < numverts; i++)
+       {
+       VECADD(verts[i].tx, verts[i].txold, verts[i].tv);
+}
+               
+               // update cloth bvh
+       bvh_update_from_cloth(clmd, 1);  // 0 means STATIC, 1 means MOVING 
+               
+               
+               // free collision list
+       if(clmd->coll_parms->collision_list)
+       {
+       LinkNode *search = clmd->coll_parms->collision_list;
+       while(search)
+       {
+       CollPair *coll_pair = search->link;
+                                                       
+       MEM_freeN(coll_pair);
+       search = search->next;
+}
+       BLI_linklist_free(clmd->coll_parms->collision_list,NULL);
+                       
+       clmd->coll_parms->collision_list = NULL;
+}
+               
+               // printf("ic: %d\n", ic);
+       rounds++;
+}
+       while(result && (CLOTH_MAX_THRESHOLD>rounds));
+       
+       ////////////////////////////////////////////////////////////
+       // update positions + velocities
+       ////////////////////////////////////////////////////////////
+       
+       // verts come from clmd
+       for(i = 0; i < numverts; i++)
+       {
+       VECADD(verts[i].tx, verts[i].txold, verts[i].tv);
+}
+       ////////////////////////////////////////////////////////////
+       */
+       
+       return MIN2(ret, 1);
+}
diff --git a/source/blender/blenkernel/intern/implicit.c b/source/blender/blenkernel/intern/implicit.c
new file mode 100644 (file)
index 0000000..f56b94c
--- /dev/null
@@ -0,0 +1,1571 @@
+/*  implicit.c      
+* 
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) Blender Foundation
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): none yet.
+*
+* ***** END GPL/BL DUAL 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 "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 "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 int first = 1;
+       if(first) {
+               QueryPerformanceFrequency(&ifreq);
+               first = 0;
+       }
+       QueryPerformanceCounter(&_itstart);
+}
+void itend(void)
+{
+       QueryPerformanceCounter(&_itend);
+}
+double itval()
+{
+       return ((double)_itend.QuadPart -
+                       (double)_itstart.QuadPart)/((double)ifreq.QuadPart);
+}
+#else
+#include <sys/time.h>
+// intrinsics need better compile flag checking
+// #include <xmmintrin.h>
+// #include <pmmintrin.h>
+// #include <pthread.h>
+
+                        static struct timeval _itstart, _itend;
+        static struct timezone itz;
+        void itstart(void)
+{
+       gettimeofday(&_itstart, &itz);
+}
+void itend(void)
+{
+       gettimeofday(&_itend,&itz);
+}
+double itval()
+{
+       double t1, t2;
+       t1 =  (double)_itstart.tv_sec + (double)_itstart.tv_usec/(1000*1000);
+       t2 =  (double)_itend.tv_sec + (double)_itend.tv_usec/(1000*1000);
+       return t2-t1;
+}
+#endif
+/*
+#define C99
+#ifdef C99
+#defineDO_INLINE inline 
+#else 
+#defineDO_INLINE static 
+#endif
+*/
+struct Cloth;
+
+//////////////////////////////////////////
+/* fast vector / matrix library, enhancements are welcome :) -dg */
+/////////////////////////////////////////
+
+/* DEFINITIONS */
+typedef float lfVector[3];
+typedef struct fmatrix3x3 {
+       float m[3][3]; /* 4x4 matrix */
+       unsigned int c,r; /* column and row number */
+       int pinned; /* is this vertex allowed to move? */
+       float n1,n2,n3; /* three normal vectors for collision constrains */
+       unsigned int vcount; /* vertex count */
+       unsigned int scount; /* spring count */ 
+} fmatrix3x3;
+
+///////////////////////////
+// float[3] vector
+///////////////////////////
+/* simple vector code */
+/* STATUS: verified */
+DO_INLINE void mul_fvector_S(float to[3], float from[3], float scalar)
+{
+       to[0] = from[0] * scalar;
+       to[1] = from[1] * scalar;
+       to[2] = from[2] * scalar;
+}
+/* simple cross product */
+/* STATUS: verified */
+DO_INLINE void cross_fvector(float to[3], float vectorA[3], float vectorB[3])
+{
+       to[0] = vectorA[1] * vectorB[2] - vectorA[2] * vectorB[1];
+       to[1] = vectorA[2] * vectorB[0] - vectorA[0] * vectorB[2];
+       to[2] = vectorA[0] * vectorB[1] - vectorA[1] * vectorB[0];
+}
+/* simple v^T * v product ("outer product") */
+/* STATUS: HAS TO BE verified (*should* work) */
+DO_INLINE void mul_fvectorT_fvector(float to[3][3], float vectorA[3], float vectorB[3])
+{
+       mul_fvector_S(to[0], vectorB, vectorA[0]);
+       mul_fvector_S(to[1], vectorB, vectorA[1]);
+       mul_fvector_S(to[2], vectorB, vectorA[2]);
+}
+/* simple v^T * v product with scalar ("outer product") */
+/* STATUS: HAS TO BE verified (*should* work) */
+DO_INLINE void mul_fvectorT_fvectorS(float to[3][3], float vectorA[3], float vectorB[3], float aS)
+{
+       mul_fvector_S(to[0], vectorB, vectorA[0]* aS);
+       mul_fvector_S(to[1], vectorB, vectorA[1]* aS);
+       mul_fvector_S(to[2], vectorB, vectorA[2]* aS);
+}
+
+/* printf vector[3] on console: for debug output */
+void print_fvector(float m3[3])
+{
+       printf("%f\n%f\n%f\n\n",m3[0],m3[1],m3[2]);
+}
+
+///////////////////////////
+// long float vector float (*)[3]
+///////////////////////////
+/* print long vector on console: for debug output */
+DO_INLINE void print_lfvector(float (*fLongVector)[3], unsigned int verts)
+{
+       unsigned int i = 0;
+       for(i = 0; i < verts; i++)
+       {
+               print_fvector(fLongVector[i]);
+       }
+}
+/* create long vector */
+DO_INLINE lfVector *create_lfvector(unsigned int verts)
+{
+       // TODO: check if memory allocation was successfull */
+       return  (lfVector *)MEM_callocN (verts * sizeof(lfVector), "cloth_implicit_alloc_vector");
+       // return (lfVector *)cloth_aligned_malloc(&MEMORY_BASE, verts * sizeof(lfVector));
+}
+/* delete long vector */
+DO_INLINE void del_lfvector(float (*fLongVector)[3])
+{
+       if (fLongVector != NULL)
+       {
+               MEM_freeN (fLongVector);
+               // cloth_aligned_free(&MEMORY_BASE, fLongVector);
+       }
+}
+/* copy long vector */
+DO_INLINE void cp_lfvector(float (*to)[3], float (*from)[3], unsigned int verts)
+{
+       memcpy(to, from, verts * sizeof(lfVector));
+}
+/* init long vector with float[3] */
+DO_INLINE void init_lfvector(float (*fLongVector)[3], float vector[3], unsigned int verts)
+{
+       unsigned int i = 0;
+       for(i = 0; i < verts; i++)
+       {
+               VECCOPY(fLongVector[i], vector);
+       }
+}
+/* zero long vector with float[3] */
+DO_INLINE void zero_lfvector(float (*to)[3], unsigned int verts)
+{
+       memset(to, 0.0f, verts * sizeof(lfVector));
+}
+/* multiply long vector with scalar*/
+DO_INLINE void mul_lfvectorS(float (*to)[3], float (*fLongVector)[3], float scalar, unsigned int verts)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < verts; i++)
+       {
+               mul_fvector_S(to[i], fLongVector[i], scalar);
+       }
+}
+/* multiply long vector with scalar*/
+/* A -= B * float */
+DO_INLINE void submul_lfvectorS(float (*to)[3], float (*fLongVector)[3], float scalar, unsigned int verts)
+{
+       unsigned int i = 0;
+       for(i = 0; i < verts; i++)
+       {
+               VECSUBMUL(to[i], fLongVector[i], scalar);
+       }
+}
+/* dot product for big vector */
+DO_INLINE float dot_lfvector(float (*fLongVectorA)[3], float (*fLongVectorB)[3], unsigned int verts)
+{
+       unsigned int i = 0;
+       float temp = 0.0;
+// schedule(guided, 2)
+#pragma omp parallel for reduction(+: temp)
+       for(i = 0; i < verts; i++)
+       {
+               temp += INPR(fLongVectorA[i], fLongVectorB[i]);
+       }
+       return temp;
+}
+/* A = B + C  --> for big vector */
+DO_INLINE void add_lfvector_lfvector(float (*to)[3], float (*fLongVectorA)[3], float (*fLongVectorB)[3], unsigned int verts)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < verts; i++)
+       {
+               VECADD(to[i], fLongVectorA[i], fLongVectorB[i]);
+       }
+
+}
+/* A = B + C * float --> for big vector */
+DO_INLINE void add_lfvector_lfvectorS(float (*to)[3], float (*fLongVectorA)[3], float (*fLongVectorB)[3], float bS, unsigned int verts)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < verts; i++)
+       {
+               VECADDS(to[i], fLongVectorA[i], fLongVectorB[i], bS);
+
+       }
+}
+/* A = B * float + C * float --> for big vector */
+DO_INLINE void add_lfvectorS_lfvectorS(float (*to)[3], float (*fLongVectorA)[3], float aS, float (*fLongVectorB)[3], float bS, unsigned int verts)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < verts; i++)
+       {
+               VECADDSS(to[i], fLongVectorA[i], aS, fLongVectorB[i], bS);
+       }
+}
+/* A = B - C * float --> for big vector */
+DO_INLINE void sub_lfvector_lfvectorS(float (*to)[3], float (*fLongVectorA)[3], float (*fLongVectorB)[3], float bS, unsigned int verts)
+{
+       unsigned int i = 0;
+       for(i = 0; i < verts; i++)
+       {
+               VECSUBS(to[i], fLongVectorA[i], fLongVectorB[i], bS);
+       }
+
+}
+/* A = B - C --> for big vector */
+DO_INLINE void sub_lfvector_lfvector(float (*to)[3], float (*fLongVectorA)[3], float (*fLongVectorB)[3], unsigned int verts)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < verts; i++)
+       {
+               VECSUB(to[i], fLongVectorA[i], fLongVectorB[i]);
+       }
+
+}
+///////////////////////////
+// 4x4 matrix
+///////////////////////////
+/* printf 4x4 matrix on console: for debug output */
+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]);
+}
+/* copy 4x4 matrix */
+DO_INLINE void cp_fmatrix(float to[3][3], float from[3][3])
+{
+       // memcpy(to, from, sizeof (float) * 9);
+       VECCOPY(to[0], from[0]);
+       VECCOPY(to[1], from[1]);
+       VECCOPY(to[2], from[2]);
+}
+/* calculate determinant of 4x4 matrix */
+DO_INLINE float det_fmatrix(float m[3][3])
+{
+       return  m[0][0]*m[1][1]*m[2][2] + m[1][0]*m[2][1]*m[0][2] + m[0][1]*m[1][2]*m[2][0] 
+                       -m[0][0]*m[1][2]*m[2][1] - m[0][1]*m[1][0]*m[2][2] - m[2][0]*m[1][1]*m[0][2];
+}
+DO_INLINE void inverse_fmatrix(float to[3][3], float from[3][3])
+{
+       unsigned int i, j;
+       float d;
+
+       if((d=det_fmatrix(from))==0)
+       {
+               printf("can't build inverse");
+               exit(0);
+       }
+       for(i=0;i<3;i++) 
+       {
+               for(j=0;j<3;j++) 
+               {
+                       int i1=(i+1)%3;
+                       int i2=(i+2)%3;
+                       int j1=(j+1)%3;
+                       int j2=(j+2)%3;
+                       // reverse indexs i&j to take transpose
+                       to[j][i] = (from[i1][j1]*from[i2][j2]-from[i1][j2]*from[i2][j1])/d;
+                       /*
+                       if(i==j)
+                       to[i][j] = 1.0f / from[i][j];
+                       else
+                       to[i][j] = 0;
+                       */
+               }
+       }
+
+}
+
+/* 4x4 matrix multiplied by a scalar */
+/* STATUS: verified */
+DO_INLINE void mul_fmatrix_S(float matrix[3][3], float scalar)
+{
+       mul_fvector_S(matrix[0], matrix[0],scalar);
+       mul_fvector_S(matrix[1], matrix[1],scalar);
+       mul_fvector_S(matrix[2], matrix[2],scalar);
+}
+
+/* a vector multiplied by a 4x4 matrix */
+/* STATUS: verified */
+DO_INLINE void mul_fvector_fmatrix(float *to, float *from, float matrix[3][3])
+{
+       to[0] = matrix[0][0]*from[0] + matrix[1][0]*from[1] + matrix[2][0]*from[2];
+       to[1] = matrix[0][1]*from[0] + matrix[1][1]*from[1] + matrix[2][1]*from[2];
+       to[2] = matrix[0][2]*from[0] + matrix[1][2]*from[1] + matrix[2][2]*from[2];
+}
+
+/* 4x4 matrix multiplied by a vector */
+/* STATUS: verified */
+DO_INLINE void mul_fmatrix_fvector(float *to, float matrix[3][3], float *from)
+{
+       to[0] = INPR(matrix[0],from);
+       to[1] = INPR(matrix[1],from);
+       to[2] = INPR(matrix[2],from);
+}
+/* 4x4 matrix multiplied by a 4x4 matrix */
+/* STATUS: verified */
+DO_INLINE void mul_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       mul_fvector_fmatrix(to[0], matrixA[0],matrixB);
+       mul_fvector_fmatrix(to[1], matrixA[1],matrixB);
+       mul_fvector_fmatrix(to[2], matrixA[2],matrixB);
+}
+/* 4x4 matrix addition with 4x4 matrix */
+DO_INLINE void add_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       VECADD(to[0], matrixA[0], matrixB[0]);
+       VECADD(to[1], matrixA[1], matrixB[1]);
+       VECADD(to[2], matrixA[2], matrixB[2]);
+}
+/* 4x4 matrix add-addition with 4x4 matrix */
+DO_INLINE void addadd_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       VECADDADD(to[0], matrixA[0], matrixB[0]);
+       VECADDADD(to[1], matrixA[1], matrixB[1]);
+       VECADDADD(to[2], matrixA[2], matrixB[2]);
+}
+/* 4x4 matrix sub-addition with 4x4 matrix */
+DO_INLINE void addsub_fmatrixS_fmatrixS(float to[3][3], float matrixA[3][3], float aS, float matrixB[3][3], float bS)
+{
+       VECADDSUBSS(to[0], matrixA[0], aS, matrixB[0], bS);
+       VECADDSUBSS(to[1], matrixA[1], aS, matrixB[1], bS);
+       VECADDSUBSS(to[2], matrixA[2], aS, matrixB[2], bS);
+}
+/* A -= B + C (4x4 matrix sub-addition with 4x4 matrix) */
+DO_INLINE void subadd_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       VECSUBADD(to[0], matrixA[0], matrixB[0]);
+       VECSUBADD(to[1], matrixA[1], matrixB[1]);
+       VECSUBADD(to[2], matrixA[2], matrixB[2]);
+}
+/* A -= B*x + C*y (4x4 matrix sub-addition with 4x4 matrix) */
+DO_INLINE void subadd_fmatrixS_fmatrixS(float to[3][3], float matrixA[3][3], float aS, float matrixB[3][3], float bS)
+{
+       VECSUBADDSS(to[0], matrixA[0], aS, matrixB[0], bS);
+       VECSUBADDSS(to[1], matrixA[1], aS, matrixB[1], bS);
+       VECSUBADDSS(to[2], matrixA[2], aS, matrixB[2], bS);
+}
+/* A = B - C (4x4 matrix subtraction with 4x4 matrix) */
+DO_INLINE void sub_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       VECSUB(to[0], matrixA[0], matrixB[0]);
+       VECSUB(to[1], matrixA[1], matrixB[1]);
+       VECSUB(to[2], matrixA[2], matrixB[2]);
+}
+/* A += B - C (4x4 matrix add-subtraction with 4x4 matrix) */
+DO_INLINE void addsub_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       VECADDSUB(to[0], matrixA[0], matrixB[0]);
+       VECADDSUB(to[1], matrixA[1], matrixB[1]);
+       VECADDSUB(to[2], matrixA[2], matrixB[2]);
+}
+/////////////////////////////////////////////////////////////////
+// special functions
+/////////////////////////////////////////////////////////////////
+/* a vector multiplied and added to/by a 4x4 matrix */
+DO_INLINE void muladd_fvector_fmatrix(float to[3], float from[3], float matrix[3][3])
+{
+       to[0] += matrix[0][0]*from[0] + matrix[1][0]*from[1] + matrix[2][0]*from[2];
+       to[1] += matrix[0][1]*from[0] + matrix[1][1]*from[1] + matrix[2][1]*from[2];
+       to[2] += matrix[0][2]*from[0] + matrix[1][2]*from[1] + matrix[2][2]*from[2];
+}
+/* 4x4 matrix multiplied and added  to/by a 4x4 matrix  and added to another 4x4 matrix */
+DO_INLINE void muladd_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       muladd_fvector_fmatrix(to[0], matrixA[0],matrixB);
+       muladd_fvector_fmatrix(to[1], matrixA[1],matrixB);
+       muladd_fvector_fmatrix(to[2], matrixA[2],matrixB);
+}
+/* a vector multiplied and sub'd to/by a 4x4 matrix */
+DO_INLINE void mulsub_fvector_fmatrix(float to[3], float from[3], float matrix[3][3])
+{
+       to[0] -= matrix[0][0]*from[0] + matrix[1][0]*from[1] + matrix[2][0]*from[2];
+       to[1] -= matrix[0][1]*from[0] + matrix[1][1]*from[1] + matrix[2][1]*from[2];
+       to[2] -= matrix[0][2]*from[0] + matrix[1][2]*from[1] + matrix[2][2]*from[2];
+}
+/* 4x4 matrix multiplied and sub'd  to/by a 4x4 matrix  and added to another 4x4 matrix */
+DO_INLINE void mulsub_fmatrix_fmatrix(float to[3][3], float matrixA[3][3], float matrixB[3][3])
+{
+       mulsub_fvector_fmatrix(to[0], matrixA[0],matrixB);
+       mulsub_fvector_fmatrix(to[1], matrixA[1],matrixB);
+       mulsub_fvector_fmatrix(to[2], matrixA[2],matrixB);
+}
+/* 4x4 matrix multiplied+added by a vector */
+/* STATUS: verified */
+DO_INLINE void muladd_fmatrix_fvector(float to[3], float matrix[3][3], float from[3])
+{
+       to[0] += INPR(matrix[0],from);
+       to[1] += INPR(matrix[1],from);
+       to[2] += INPR(matrix[2],from);  
+}
+/* 4x4 matrix multiplied+sub'ed by a vector */
+DO_INLINE void mulsub_fmatrix_fvector(float to[3], float matrix[3][3], float from[3])
+{
+       to[0] -= INPR(matrix[0],from);
+       to[1] -= INPR(matrix[1],from);
+       to[2] -= INPR(matrix[2],from);
+}
+/////////////////////////////////////////////////////////////////
+
+///////////////////////////
+// SPARSE SYMMETRIC big matrix with 4x4 matrix entries
+///////////////////////////
+/* printf a big matrix on console: for debug output */
+void print_bfmatrix(fmatrix3x3 *m3)
+{
+       unsigned int i = 0;
+
+       for(i = 0; i < m3[0].vcount + m3[0].scount; i++)
+       {
+               print_fmatrix(m3[i].m);
+       }
+}
+/* create big matrix */
+DO_INLINE fmatrix3x3 *create_bfmatrix(unsigned int verts, unsigned int springs)
+{
+       // TODO: check if memory allocation was successfull */
+       fmatrix3x3 *temp = (fmatrix3x3 *)MEM_callocN (sizeof (fmatrix3x3) * (verts + springs), "cloth_implicit_alloc_matrix");
+       temp[0].vcount = verts;
+       temp[0].scount = springs;
+       return temp;
+}
+/* delete big matrix */
+DO_INLINE void del_bfmatrix(fmatrix3x3 *matrix)
+{
+       if (matrix != NULL)
+       {
+               MEM_freeN (matrix);
+       }
+}
+/* copy big matrix */
+DO_INLINE void cp_bfmatrix(fmatrix3x3 *to, fmatrix3x3 *from)
+{      
+       // TODO bounds checking 
+       memcpy(to, from, sizeof(fmatrix3x3) * (from[0].vcount+from[0].scount) );
+}
+/* init the diagonal of big matrix */
+// slow in parallel
+DO_INLINE void initdiag_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
+{
+       unsigned int i,j;
+       float tmatrix[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
+
+       for(i = 0; i < matrix[0].vcount; i++)
+       {               
+               cp_fmatrix(matrix[i].m, m3); 
+       }
+       for(j = matrix[0].vcount; j < matrix[0].vcount+matrix[0].scount; j++)
+       {
+               cp_fmatrix(matrix[j].m, tmatrix); 
+       }
+}
+/* init big matrix */
+DO_INLINE void init_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
+{
+       unsigned int i;
+
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               cp_fmatrix(matrix[i].m, m3); 
+       }
+}
+/* multiply big matrix with scalar*/
+DO_INLINE void mul_bfmatrix_S(fmatrix3x3 *matrix, float scalar)
+{
+       unsigned int i = 0;
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               mul_fmatrix_S(matrix[i].m, scalar);
+       }
+}
+
+/* SPARSE SYMMETRIC multiply big matrix with long vector*/
+/* STATUS: verified */
+DO_INLINE void mul_bfmatrix_lfvector( float (*to)[3], fmatrix3x3 *from, lfVector *fLongVector)
+{
+       unsigned int i = 0;
+       lfVector *temp = create_lfvector(from[0].vcount);
+       
+       zero_lfvector(to, from[0].vcount);
+
+#pragma omp parallel sections private(i)
+       {
+#pragma omp section
+               {
+                       for(i = from[0].vcount; i < from[0].vcount+from[0].scount; i++)
+                       {
+                               muladd_fmatrix_fvector(to[from[i].c], from[i].m, fLongVector[from[i].r]);
+                       }
+               }       
+#pragma omp section
+               {
+                       for(i = 0; i < from[0].vcount+from[0].scount; i++)
+                       {
+                               muladd_fmatrix_fvector(temp[from[i].r], from[i].m, fLongVector[from[i].c]);
+                       }
+               }
+       }
+       add_lfvector_lfvector(to, to, temp, from[0].vcount);
+       
+       del_lfvector(temp);
+       
+       
+}
+
+/* SPARSE SYMMETRIC multiply big matrix with long vector (for diagonal preconditioner) */
+/* STATUS: verified */
+DO_INLINE void mul_prevfmatrix_lfvector( float (*to)[3], fmatrix3x3 *from, lfVector *fLongVector)
+{
+       unsigned int i = 0;
+       
+       for(i = 0; i < from[0].vcount; i++)
+       {
+               mul_fmatrix_fvector(to[from[i].r], from[i].m, fLongVector[from[i].c]);
+       }
+}
+
+/* SPARSE SYMMETRIC add big matrix with big matrix: A = B + C*/
+DO_INLINE void add_bfmatrix_bfmatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               add_fmatrix_fmatrix(to[i].m, from[i].m, matrix[i].m);   
+       }
+
+}
+/* SPARSE SYMMETRIC add big matrix with big matrix: A += B + C */
+DO_INLINE void addadd_bfmatrix_bfmatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               addadd_fmatrix_fmatrix(to[i].m, from[i].m, matrix[i].m);        
+       }
+
+}
+/* SPARSE SYMMETRIC subadd big matrix with big matrix: A -= B + C */
+DO_INLINE void subadd_bfmatrix_bfmatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               subadd_fmatrix_fmatrix(to[i].m, from[i].m, matrix[i].m);        
+       }
+
+}
+/*  A = B - C (SPARSE SYMMETRIC sub big matrix with big matrix) */
+DO_INLINE void sub_bfmatrix_bfmatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               sub_fmatrix_fmatrix(to[i].m, from[i].m, matrix[i].m);   
+       }
+
+}
+/* SPARSE SYMMETRIC sub big matrix with big matrix S (special constraint matrix with limited entries) */
+DO_INLINE void sub_bfmatrix_Smatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount; i++)
+       {
+               sub_fmatrix_fmatrix(to[matrix[i].c].m, from[matrix[i].c].m, matrix[i].m);       
+       }
+
+}
+/* A += B - C (SPARSE SYMMETRIC addsub big matrix with big matrix) */
+DO_INLINE void addsub_bfmatrix_bfmatrix( fmatrix3x3 *to, fmatrix3x3 *from,  fmatrix3x3 *matrix)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               addsub_fmatrix_fmatrix(to[i].m, from[i].m, matrix[i].m);        
+       }
+
+}
+/* SPARSE SYMMETRIC sub big matrix with big matrix*/
+/* A -= B * float + C * float --> for big matrix */
+/* VERIFIED */
+DO_INLINE void subadd_bfmatrixS_bfmatrixS( fmatrix3x3 *to, fmatrix3x3 *from, float aS,  fmatrix3x3 *matrix, float bS)
+{
+       unsigned int i = 0;
+
+       /* process diagonal elements */
+       for(i = 0; i < matrix[0].vcount+matrix[0].scount; i++)
+       {
+               subadd_fmatrixS_fmatrixS(to[i].m, from[i].m, aS, matrix[i].m, bS);      
+       }
+
+}
+
+///////////////////////////////////////////////////////////////////
+// simulator start
+///////////////////////////////////////////////////////////////////
+static float I[3][3] = {{1,0,0},{0,1,0},{0,0,1}};
+static float ZERO[3][3] = {{0,0,0}, {0,0,0}, {0,0,0}};
+typedef struct Implicit_Data 
+{
+       lfVector *X, *V, *Xnew, *Vnew, *olddV, *F, *B, *dV, *z;
+       fmatrix3x3 *A, *dFdV, *dFdX, *S, *P, *Pinv, *bigI; 
+} Implicit_Data;
+
+int implicit_init (Object *ob, ClothModifierData *clmd)
+{
+       unsigned int i = 0;
+       unsigned int pinned = 0;
+       Cloth *cloth = NULL;
+       ClothVertex *verts = NULL;
+       ClothSpring *spring = NULL;
+       Implicit_Data *id = NULL;
+       LinkNode *search = NULL;
+       
+       if(G.rt > 0)
+               printf("implicit_init\n");
+
+       // init memory guard
+       // MEMORY_BASE.first = MEMORY_BASE.last = NULL;
+
+       cloth = (Cloth *)clmd->clothObject;
+       verts = cloth->verts;
+
+       // create implicit base
+       id = (Implicit_Data *)MEM_callocN (sizeof(Implicit_Data), "implicit vecmat");
+       cloth->implicit = id;
+
+       /* process diagonal elements */         
+       id->A = create_bfmatrix(cloth->numverts, cloth->numsprings);
+       id->dFdV = create_bfmatrix(cloth->numverts, cloth->numsprings);
+       id->dFdX = create_bfmatrix(cloth->numverts, cloth->numsprings);
+       id->S = create_bfmatrix(cloth->numverts, 0);
+       id->Pinv = create_bfmatrix(cloth->numverts, cloth->numsprings);
+       id->P = create_bfmatrix(cloth->numverts, cloth->numsprings);
+       id->bigI = create_bfmatrix(cloth->numverts, cloth->numsprings); // TODO 0 springs
+       id->X = create_lfvector(cloth->numverts);
+       id->Xnew = create_lfvector(cloth->numverts);
+       id->V = create_lfvector(cloth->numverts);
+       id->Vnew = create_lfvector(cloth->numverts);
+       id->olddV = create_lfvector(cloth->numverts);
+       zero_lfvector(id->olddV, cloth->numverts);
+       id->F = create_lfvector(cloth->numverts);
+       id->B = create_lfvector(cloth->numverts);
+       id->dV = create_lfvector(cloth->numverts);
+       id->z = create_lfvector(cloth->numverts);
+       
+       for(i=0;i<cloth->numverts;i++) 
+       {
+               id->A[i].r = id->A[i].c = id->dFdV[i].r = id->dFdV[i].c = id->dFdX[i].r = id->dFdX[i].c = id->P[i].c = id->P[i].r = id->Pinv[i].c = id->Pinv[i].r = id->bigI[i].c = id->bigI[i].r = i;
+
+               if(verts [i].goal >= SOFTGOALSNAP)
+               {
+                       id->S[pinned].pinned = 1;
+                       id->S[pinned].c = id->S[pinned].r = i;
+                       pinned++;
+               }
+       }
+
+       // S is special and needs specific vcount and scount
+       id->S[0].vcount = pinned; id->S[0].scount = 0;
+
+       // init springs */
+       search = cloth->springs;
+       for(i=0;i<cloth->numsprings;i++) 
+       {
+               spring = search->link;
+               
+               // dFdV_start[i].r = big_I[i].r = big_zero[i].r = 
+               id->A[i+cloth->numverts].r = id->dFdV[i+cloth->numverts].r = id->dFdX[i+cloth->numverts].r = 
+                               id->P[i+cloth->numverts].r = id->Pinv[i+cloth->numverts].r = id->bigI[i+cloth->numverts].r = spring->ij;
+
+               // dFdV_start[i].c = big_I[i].c = big_zero[i].c = 
+               id->A[i+cloth->numverts].c = id->dFdV[i+cloth->numverts].c = id->dFdX[i+cloth->numverts].c = 
+                               id->P[i+cloth->numverts].c = id->Pinv[i+cloth->numverts].c = id->bigI[i+cloth->numverts].c = spring->kl;
+
+               spring->matrix_index = i + cloth->numverts;
+               
+               search = search->next;
+       }
+
+       for(i = 0; i < cloth->numverts; i++)
+       {               
+               VECCOPY(id->X[i], verts[i].x);
+       }
+
+       return 1;
+}
+int    implicit_free (ClothModifierData *clmd)
+{
+       Implicit_Data *id;
+       Cloth *cloth;
+       cloth = (Cloth *)clmd->clothObject;
+
+       if(cloth)
+       {
+               id = cloth->implicit;
+
+               if(id)
+               {
+                       del_bfmatrix(id->A);
+                       del_bfmatrix(id->dFdV);
+                       del_bfmatrix(id->dFdX);
+                       del_bfmatrix(id->S);
+                       del_bfmatrix(id->P);
+                       del_bfmatrix(id->Pinv);
+                       del_bfmatrix(id->bigI);
+
+                       del_lfvector(id->X);
+                       del_lfvector(id->Xnew);
+                       del_lfvector(id->V);
+                       del_lfvector(id->Vnew);
+                       del_lfvector(id->olddV);
+                       del_lfvector(id->F);
+                       del_lfvector(id->B);
+                       del_lfvector(id->dV);
+                       del_lfvector(id->z);
+
+                       MEM_freeN(id);
+               }
+       }
+
+       return 1;
+}
+
+DO_INLINE float fb(float length, float L)
+{
+       float x = length/L;
+       return (-11.541f*pow(x,4)+34.193f*pow(x,3)-39.083f*pow(x,2)+23.116f*x-9.713f);
+}
+
+DO_INLINE float fbderiv(float length, float L)
+{
+       float x = length/L;
+
+       return (-46.164f*pow(x,3)+102.579f*pow(x,2)-78.166f*x+23.116f);
+}
+
+DO_INLINE float fbstar(float length, float L, float kb, float cb)
+{
+       float tempfb = kb * fb(length, L);
+
+       float fbstar = cb * (length - L);
+
+       if(tempfb < fbstar)
+               return fbstar;
+       else
+               return tempfb;          
+}
+
+// function to calculae bending spring force (taken from Choi & Co)
+DO_INLINE float fbstar_jacobi(float length, float L, float kb, float cb)
+{
+       float tempfb = kb * fb(length, L);
+       float fbstar = cb * (length - L);
+
+       if(tempfb < fbstar)
+       {               
+               return cb;
+       }
+       else
+       {
+               return kb * fbderiv(length, L); 
+       }       
+}
+
+DO_INLINE void filter(lfVector *V, fmatrix3x3 *S)
+{
+       unsigned int i=0;
+
+       for(i=0;i<S[0].vcount;i++)
+       {
+               mul_fvector_fmatrix(V[S[i].r], V[S[i].r], S[i].m);
+       }
+}
+
+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;
+       float conjgrad_epsilon=0.0001f, conjgrad_lasterror=0;
+       lfVector *q, *d, *tmp, *r; 
+       float s, starget, a, s_prev;
+       unsigned int numverts = lA[0].vcount;
+       q = create_lfvector(numverts);
+       d = create_lfvector(numverts);
+       tmp = create_lfvector(numverts);
+       r = create_lfvector(numverts);
+
+       // zero_lfvector(ldV, CLOTHPARTICLES);
+       filter(ldV, S);
+
+       add_lfvector_lfvector(ldV, ldV, z, numverts);
+
+       // r = B - Mul(tmp,A,X);    // just use B if X known to be zero
+       cp_lfvector(r, lB, numverts);
+       mul_bfmatrix_lfvector(tmp, lA, ldV);
+       sub_lfvector_lfvector(r, r, tmp, numverts);
+
+       filter(r,S);
+
+       cp_lfvector(d, r, numverts);
+
+       s = dot_lfvector(r, r, numverts);
+       starget = s * sqrt(conjgrad_epsilon);
+
+       while((s>starget && conjgrad_loopcount < conjgrad_looplimit))
+       {       
+               // Mul(q,A,d); // q = A*d;
+               mul_bfmatrix_lfvector(q, lA, d);
+
+               filter(q,S);
+
+               a = s/dot_lfvector(d, q, numverts);
+
+               // X = X + d*a;
+               add_lfvector_lfvectorS(ldV, ldV, d, a, numverts);
+
+               // r = r - q*a;
+               sub_lfvector_lfvectorS(r, r, q, a, numverts);
+
+               s_prev = s;
+               s = dot_lfvector(r, r, numverts);
+
+               //d = r+d*(s/s_prev);
+               add_lfvector_lfvectorS(d, r, d, (s/s_prev), numverts);
+
+               filter(d,S);
+
+               conjgrad_loopcount++;
+       }
+       conjgrad_lasterror = s;
+
+       del_lfvector(q);
+       del_lfvector(d);
+       del_lfvector(tmp);
+       del_lfvector(r);
+       // printf("W/O conjgrad_loopcount: %d\n", conjgrad_loopcount);
+
+       return conjgrad_loopcount<conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
+}
+
+// block diagonalizer
+DO_INLINE void BuildPPinv(fmatrix3x3 *lA, fmatrix3x3 *P, fmatrix3x3 *Pinv)
+{
+       unsigned int i = 0;
+       
+       // Take only the diagonal blocks of A
+// #pragma omp parallel for private(i)
+       for(i = 0; i<lA[0].vcount; i++)
+       {
+               // block diagonalizer
+               cp_fmatrix(P[i].m, lA[i].m);
+               inverse_fmatrix(Pinv[i].m, P[i].m);
+               
+       }
+}
+
+// version 1.3
+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;
+       float conjgrad_epsilon=0.0001; // 0.2 is dt for steps=5
+       lfVector *r = create_lfvector(numverts);
+       lfVector *p = create_lfvector(numverts);
+       lfVector *s = create_lfvector(numverts);
+       lfVector *h = create_lfvector(numverts);
+       
+       BuildPPinv(lA, P, Pinv);
+       
+       filter(dv, S);
+       add_lfvector_lfvector(dv, dv, z, numverts);
+       
+       mul_bfmatrix_lfvector(r, lA, dv);
+       sub_lfvector_lfvector(r, lB, r, numverts);
+       filter(r, S);
+       
+       mul_prevfmatrix_lfvector(p, Pinv, r);
+       filter(p, S);
+       
+       deltaNew = dot_lfvector(r, p, numverts);
+       
+       delta0 = deltaNew * sqrt(conjgrad_epsilon);
+       
+       // itstart();
+       
+       while ((deltaNew > delta0) && (iterations < conjgrad_looplimit))
+       {
+               iterations++;
+               
+               mul_bfmatrix_lfvector(s, lA, p);
+               filter(s, S);
+               
+               alpha = deltaNew / dot_lfvector(p, s, numverts);
+               
+               add_lfvector_lfvectorS(dv, dv, p, alpha, numverts);
+               
+               add_lfvector_lfvectorS(r, r, s, -alpha, numverts);
+               
+               mul_prevfmatrix_lfvector(h, Pinv, r);
+               filter(h, S);
+               
+               deltaOld = deltaNew;
+               
+               deltaNew = dot_lfvector(r, h, numverts);
+               
+               add_lfvector_lfvectorS(p, h, p, deltaNew / deltaOld, numverts);
+               
+               filter(p, S);
+               
+       }
+       
+       // itend();
+       // printf("cg_filtered_pre time: %f\n", (float)itval());
+       
+       del_lfvector(h);
+       del_lfvector(s);
+       del_lfvector(p);
+       del_lfvector(r);
+       
+       // printf("iterations: %d\n", iterations);
+       
+       return iterations<conjgrad_looplimit;
+}
+
+// outer product is NOT cross product!!!
+DO_INLINE void dfdx_spring_type1(float to[3][3], float dir[3],float length,float L,float k)
+{
+       // dir is unit length direction, rest is spring's restlength, k is spring constant.
+       // return  (outerprod(dir,dir)*k + (I - outerprod(dir,dir))*(k - ((k*L)/length)));
+       float temp[3][3];
+       mul_fvectorT_fvector(temp, dir, dir);
+       sub_fmatrix_fmatrix(to, I, temp);
+       mul_fmatrix_S(to, k* (1.0f-(L/length)));
+       mul_fmatrix_S(temp, k);
+       add_fmatrix_fmatrix(to, temp, to);
+}
+
+DO_INLINE void dfdx_spring_type2(float to[3][3], float dir[3],float length,float L,float k, float cb)
+{
+       // return  outerprod(dir,dir)*fbstar_jacobi(length, L, k, cb);
+       mul_fvectorT_fvectorS(to, dir, dir, fbstar_jacobi(length, L, k, cb));
+}
+
+DO_INLINE void dfdv_damp(float to[3][3], float dir[3], float damping)
+{
+       // derivative of force wrt velocity.  
+       // return outerprod(dir,dir) * damping;
+       mul_fvectorT_fvectorS(to, dir, dir, damping);
+}
+
+DO_INLINE void dfdx_spring(float to[3][3],  float dir[3],float length,float L,float k)
+{
+       // dir is unit length direction, rest is spring's restlength, k is spring constant.
+       //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))); 
+       sub_fmatrix_fmatrix(to, to, I);
+       mul_fmatrix_S(to, -k);
+}
+
+DO_INLINE void dfdx_damp(float to[3][3],  float dir[3],float length,const float vel[3],float rest,float damping)
+{
+       // inner spring damping   vel is the relative velocity  of the endpoints.  
+       //      return (I-outerprod(dir,dir)) * (-damping * -(dot(dir,vel)/Max(length,rest)));
+       mul_fvectorT_fvector(to, dir, dir);
+       sub_fmatrix_fmatrix(to, I, to);
+       mul_fmatrix_S(to,  (-damping * -(INPR(dir,vel)/MAX2(length,rest)))); 
+
+}
+
+DO_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s, lfVector *lF, lfVector *X, lfVector *V, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX)
+{
+       float extent[3];
+       float length = 0;
+       float dir[3] = {0,0,0};
+       float vel[3];
+       float k = 0.0f;
+       float L = s->restlen;
+       float cb = clmd->sim_parms->structural;
+
+       float nullf[3] = {0,0,0};
+       float stretch_force[3] = {0,0,0};
+       float bending_force[3] = {0,0,0};
+       float damping_force[3] = {0,0,0};
+       float nulldfdx[3][3]={ {0,0,0}, {0,0,0}, {0,0,0}};
+       
+       float scaling = 0.0;
+       
+       VECCOPY(s->f, nullf);
+       cp_fmatrix(s->dfdx, nulldfdx);
+       cp_fmatrix(s->dfdv, nulldfdx);
+
+       // calculate elonglation
+       VECSUB(extent, X[s->kl], X[s->ij]);
+       VECSUB(vel, V[s->kl], V[s->ij]);
+       length = sqrt(INPR(extent, extent));
+       
+       s->flags &= ~CLOTH_SPRING_FLAG_NEEDED;
+       
+       if(length > ABS(ALMOST_ZERO))
+       {
+               /*
+               if(length>L)
+               {
+               if((clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED) 
+               && ((((length-L)*100.0f/L) > clmd->sim_parms->maxspringlen))) // cut spring!
+               {
+               s->flags |= CSPRING_FLAG_DEACTIVATE;
+               return;
+       }
+       } 
+               */
+               mul_fvector_S(dir, extent, 1.0f/length);
+       }
+       else    
+       {
+               mul_fvector_S(dir, extent, 0.0f);
+       }
+       
+       // calculate force of structural + shear springs
+       if(s->type != CLOTH_SPRING_TYPE_BENDING)
+       {
+               if(length > L) // only on elonglation
+               {
+                       s->flags |= CLOTH_SPRING_FLAG_NEEDED;
+                       
+                       k = clmd->sim_parms->structural;
+                               
+                       scaling = k + s->stiffness * ABS(clmd->sim_parms->max_struct-k);
+                       
+                       k = scaling / (clmd->coll_parms->avg_spring_len + FLT_EPSILON);
+                       
+                       // printf("scaling: %f, stiffness: %f\n", k, s->stiffness);
+                       
+                       /*
+                       if((s->ij == 109) || (s->kl == 109))
+                       {
+                       printf("length-L: %f, f: %f, len: %f, L: %f\n", length-L, (k*(length-L)), length, L);
+                       printf("kl X-x: %f, f-y: %f, f-z: %f\n", X[s->kl][0], X[s->kl][1], X[s->kl][2]);
+                       printf("ij X-x: %f, f-y: %f, f-z: %f\n\n", X[s->ij][0], X[s->ij][1], X[s->ij][2]);
+               }
+                       */
+
+                       mul_fvector_S(stretch_force, dir, (k*(length-L))); 
+
+                       VECADD(s->f, s->f, stretch_force);
+
+                       // Ascher & Boxman, p.21: Damping only during elonglation
+                       mul_fvector_S(damping_force, extent, clmd->sim_parms->Cdis * ((INPR(vel,extent)/length))); 
+                       VECADD(s->f, s->f, damping_force);
+
+                       dfdx_spring_type1(s->dfdx, dir,length,L,k);
+
+                       dfdv_damp(s->dfdv, dir,clmd->sim_parms->Cdis);
+               }
+       }
+       else // calculate force of bending springs
+       {
+               if(length < L)
+               {
+                       s->flags |= CLOTH_SPRING_FLAG_NEEDED;
+                       
+                       k = clmd->sim_parms->bending;   
+                       
+                       scaling = k + s->stiffness * ABS(clmd->sim_parms->max_bend-k);                  
+                       cb = k = scaling / (20.0*(clmd->coll_parms->avg_spring_len + FLT_EPSILON));
+                       
+                       // printf("scaling: %f, stiffness: %f\n", k, s->stiffness);
+
+                       mul_fvector_S(bending_force, dir, fbstar(length, L, k, cb));
+                       VECADD(s->f, s->f, bending_force);
+
+                       dfdx_spring_type2(s->dfdx, dir,length,L,k, cb);
+               }
+       }
+       /*
+       if((s->ij == 109) || (s->kl == 109))
+       {
+       printf("type: %d, f-x: %f, f-y: %f, f-z: %f\n", s->type, s->f[0], s->f[1], s->f[2]);
+}
+       */
+}
+
+DO_INLINE void cloth_apply_spring_force(ClothModifierData *clmd, ClothSpring *s, lfVector *lF, lfVector *X, lfVector *V, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX)
+{
+       if(s->flags & CLOTH_SPRING_FLAG_NEEDED)
+       {
+               if(s->type != CLOTH_SPRING_TYPE_BENDING)
+               {
+                       sub_fmatrix_fmatrix(dFdV[s->ij].m, dFdV[s->ij].m, s->dfdv);
+                       sub_fmatrix_fmatrix(dFdV[s->kl].m, dFdV[s->kl].m, s->dfdv);
+                       add_fmatrix_fmatrix(dFdV[s->matrix_index].m, dFdV[s->matrix_index].m, s->dfdv); 
+               }
+
+               VECADD(lF[s->ij], lF[s->ij], s->f);
+               VECSUB(lF[s->kl], lF[s->kl], s->f);
+
+               sub_fmatrix_fmatrix(dFdX[s->ij].m, dFdX[s->ij].m, s->dfdx);
+               sub_fmatrix_fmatrix(dFdX[s->kl].m, dFdX[s->kl].m, s->dfdx);
+
+               add_fmatrix_fmatrix(dFdX[s->matrix_index].m, dFdX[s->matrix_index].m, s->dfdx);
+       }       
+}
+
+DO_INLINE void calculateTriangleNormal(float to[3], lfVector *X, MFace mface)
+{
+       float v1[3], v2[3];
+
+       VECSUB(v1, X[mface.v2], X[mface.v1]);
+       VECSUB(v2, X[mface.v3], X[mface.v1]);
+       cross_fvector(to, v1, v2);
+}
+
+DO_INLINE void calculatQuadNormal(float to[3], lfVector *X, MFace mface)
+{
+       float temp = CalcNormFloat4(X[mface.v1],X[mface.v2],X[mface.v3],X[mface.v4],to);
+       mul_fvector_S(to, to, temp);
+}
+
+void calculateWeightedVertexNormal(ClothModifierData *clmd, MFace *mfaces, float to[3], int index, lfVector *X)
+{
+       float temp[3]; 
+       int i;
+       Cloth *cloth = clmd->clothObject;
+
+       for(i = 0; i < cloth->numfaces; i++)
+       {
+               // 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);
+               }
+       }
+}
+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)
+{      
+
+}
+
+void cloth_calc_force(ClothModifierData *clmd, lfVector *lF, lfVector *lX, lfVector *lV, fmatrix3x3 *dFdV, fmatrix3x3 *dFdX, ListBase *effectors, float time)
+{
+       /* Collect forces and derivatives:  F,dFdX,dFdV */
+       Cloth           *cloth          = clmd->clothObject;
+       unsigned int    i               = 0;
+       float           spring_air      = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
+       float           gravity[3];
+       float           tm2[3][3]       = {{-spring_air,0,0}, {0,-spring_air,0},{0,0,-spring_air}};
+       ClothVertex *verts = cloth->verts;
+       MFace           *mfaces         = cloth->mfaces;
+       float wind_normalized[3];
+       unsigned int numverts = cloth->numverts;
+       float auxvect[3], velgoal[3], tvect[3];
+       float kd, ks;
+       LinkNode *search = cloth->springs;
+
+
+       VECCOPY(gravity, clmd->sim_parms->gravity);
+       mul_fvector_S(gravity, gravity, 0.001f); /* scale gravity force */
+
+       /* set dFdX jacobi matrix to zero */
+       init_bfmatrix(dFdX, ZERO);
+       /* set dFdX jacobi matrix diagonal entries to -spring_air */ 
+       initdiag_bfmatrix(dFdV, tm2);
+
+       init_lfvector(lF, gravity, numverts);
+
+       submul_lfvectorS(lF, lV, spring_air, numverts);
+               
+       /* do goal stuff */
+       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) 
+       {       
+               for(i = 0; i < numverts; i++)
+               {                       
+                       if(verts [i].goal < SOFTGOALSNAP)
+                       {                       
+                               // current_position = xold + t * (newposition - xold)
+                               VECSUB(tvect, verts[i].xconst, verts[i].xold);
+                               mul_fvector_S(tvect, tvect, time);
+                               VECADD(tvect, tvect, verts[i].xold);
+
+                               VECSUB(auxvect, tvect, lX[i]);
+                               ks  = 1.0f/(1.0f- verts [i].goal*clmd->sim_parms->goalspring)-1.0f ;
+                               VECADDS(lF[i], lF[i], auxvect, -ks);
+
+                               // calulate damping forces generated by goals
+                               
+                               VECSUB(velgoal,verts[i].xold, verts[i].xconst);
+                               kd =  clmd->sim_parms->goalfrict * 0.01f; // friction force scale taken from SB
+                               VECSUBADDSS(lF[i], velgoal, kd, lV[i], kd);
+                               
+                               
+                       }
+               }       
+       }
+       
+       /* 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};
+               
+#pragma omp parallel for private (i) shared(lF)
+               for(i = 0; i < cloth->numverts; i++)
+               {
+                       float vertexnormal[3]={0,0,0};
+                       float fieldfactor = 1000.0f; // windfactor  = 250.0f; // from sb
+                       
+                       pdDoEffectors(effectors, lX[i], force, speed, (float)G.scene->r.cfra, 0.0f, PE_WIND_AS_SPEED);          
+                       
+                       // TODO apply forcefields here
+                       VECADDS(lF[i], lF[i], force, fieldfactor*0.01f);
+
+                       VECCOPY(wind_normalized, speed);
+                       Normalize(wind_normalized);
+                       
+                       calculateWeightedVertexNormal(clmd, mfaces, vertexnormal, i, lX);
+                       VECADDS(lF[i], lF[i], wind_normalized, -calculateVertexWindForce(speed, vertexnormal));
+               }
+       }
+               
+       // calculate spring forces
+       search = cloth->springs;
+       while(search)
+       {
+               // only handle active springs
+               // if(((clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED) && !(springs[i].flags & CSPRING_FLAG_DEACTIVATE))|| !(clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED)){}
+               cloth_calc_spring_force(clmd, search->link, lF, lX, lV, dFdV, dFdX);
+
+               search = search->next;
+       }
+       
+       // apply spring forces
+       search = cloth->springs;
+       while(search)
+       {
+               // only handle active springs
+               // if(((clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED) && !(springs[i].flags & CSPRING_FLAG_DEACTIVATE))|| !(clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED))  
+               cloth_apply_spring_force(clmd, search->link, lF, lX, lV, dFdV, dFdX);
+               search = search->next;
+       }
+       // 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)
+{
+       unsigned int numverts = dFdV[0].vcount;
+
+       lfVector *dFdXmV = create_lfvector(numverts);
+       initdiag_bfmatrix(A, I);
+       zero_lfvector(dV, numverts);
+
+       subadd_bfmatrixS_bfmatrixS(A, dFdV, dt, dFdX, (dt*dt));
+
+       mul_bfmatrix_lfvector(dFdXmV, dFdX, lV);
+
+       add_lfvectorS_lfvectorS(B, lF, dt, dFdXmV, (dt*dt), numverts);
+       
+       itstart();
+       
+       cg_filtered(dV, A, B, z, S); /* conjugate gradient algorithm to solve Ax=b */
+       // cg_filtered_pre(dV, A, B, z, olddV, P, Pinv, dt);
+       
+       itend();
+       // printf("cg_filtered calc time: %f\n", (float)itval());
+       
+       cp_lfvector(olddV, dV, numverts);
+
+       // advance velocities
+       add_lfvector_lfvector(Vnew, lV, dV, numverts);
+       
+
+       del_lfvector(dFdXmV);
+}
+
+int implicit_solver (Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors)
+{              
+       unsigned int i=0;
+       float step=0.0f, tf=1.0f;
+       Cloth *cloth = clmd->clothObject;
+       ClothVertex *verts = cloth->verts;
+       unsigned int numverts = cloth->numverts;
+       float dt = 1.0f / clmd->sim_parms->stepsPerFrame;
+       Implicit_Data *id = cloth->implicit;
+       int result = 0;
+       
+       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) /* do goal stuff */
+       {
+               for(i = 0; i < numverts; i++)
+               {                       
+                       // update velocities with constrained velocities from pinned verts
+                       if(verts [i].goal >= SOFTGOALSNAP)
+                       {                       
+                               VECSUB(id->V[i], verts[i].xconst, verts[i].xold);
+                               // VecMulf(id->V[i], 1.0 / dt);
+                       }
+               }       
+       }
+
+       while(step < tf)
+       {               
+               effectors= pdInitEffectors(ob,NULL);
+               
+               // calculate 
+               cloth_calc_force(clmd, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step );
+               
+               // printf("F -> x: %f, y: %f; z: %f\n\n", id->F[109][0], id->F[109][1], id->F[109][2]);
+                       
+               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);
+               
+               add_lfvector_lfvectorS(id->Xnew, id->X, id->Vnew, dt, numverts);
+               
+               /*
+               printf("dt: %f\n", dt);
+               printf("Xnew -> x: %f, y: %f; z: %f\n", id->Xnew[109][0], id->Xnew[109][1], id->Xnew[109][2]);
+               printf("X    -> x: %f, y: %f; z: %f\n", id->X[109][0], id->X[109][1], id->X[109][2]);
+               printf("Vnew -> x: %f, y: %f; z: %f\n\n", id->Vnew[109][0], id->Vnew[109][1], id->Vnew[109][2]);
+               */
+               
+               // clmd->coll_parms->flags &= ~CLOTH_COLLSETTINGS_FLAG_ENABLED;
+               
+               if(clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED)
+               {
+                       // collisions 
+                       // itstart();
+                       
+                       // update verts to current positions
+                       for(i = 0; i < numverts; i++)
+                       {               
+                               
+                               if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) 
+                               {                       
+                                       if(verts [i].goal >= SOFTGOALSNAP)
+                                       {                       
+                                               float tvect[3] = {.0,.0,.0};
+                                               // VECSUB(tvect, id->Xnew[i], verts[i].xold);
+                                               mul_fvector_S(tvect, id->V[i], step+dt);
+                                               VECADD(tvect, tvect, verts[i].xold);
+                                               VECCOPY(id->Xnew[i], tvect);
+                                       }
+                                               
+                               }
+                               
+                               VECCOPY(verts[i].tx, id->Xnew[i]);
+                               
+                               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);
+       
+                       // copy corrected positions back to simulation
+                       for(i = 0; i < numverts; i++)
+                       {               
+                               if(result)
+                               {
+                                       
+                                       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) 
+                                       {                       
+                                               if(verts [i].goal >= SOFTGOALSNAP)
+                                               {
+                                                       continue;
+                                               }
+                                       }
+                                       
+                                               
+                                       // VECADD(verts[i].tx, verts[i].txold, verts[i].tv);
+                                       
+                                       VECCOPY(verts[i].txold, verts[i].tx);
+                                       
+                                       VECCOPY(id->Xnew[i], verts[i].tx);
+                                       
+                                       VECCOPY(id->Vnew[i], verts[i].tv);
+                                       VecMulf(id->Vnew[i], 1.0f / dt);
+                               }
+                               else
+                               {
+                                       VECCOPY(verts[i].txold, id->Xnew[i]);
+                               }
+                       }
+                       
+                       // X = Xnew;
+                       cp_lfvector(id->X, id->Xnew, numverts);
+                       
+                       // if there were collisions, advance the velocity from v_n+1/2 to v_n+1
+                       if(result)
+                       {
+                               // V = Vnew;
+                               cp_lfvector(id->V, id->Vnew, numverts);
+                               
+                               // calculate 
+                               cloth_calc_force(clmd, id->F, id->X, id->V, id->dFdV, id->dFdX, effectors, step);       
+                               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);
+                       }
+               }
+               else
+               {
+                       // X = Xnew;
+                       cp_lfvector(id->X, id->Xnew, numverts);
+               }
+               
+               // itend();
+               // printf("collision time: %f\n", (float)itval());
+               
+               // V = Vnew;
+               cp_lfvector(id->V, id->Vnew, numverts);
+
+               step += dt;
+
+               if(effectors) pdEndEffectors(effectors);
+       }
+
+       for(i = 0; i < numverts; i++)
+       {                               
+               if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL)
+               {
+                       if(verts [i].goal < SOFTGOALSNAP)
+                       {
+                               VECCOPY(verts[i].txold, id->X[i]);
+                               VECCOPY(verts[i].x, id->X[i]);
+                               VECCOPY(verts[i].v, id->V[i]);
+                       }
+                       else
+                       {
+                               VECCOPY(verts[i].txold, verts[i].xconst);
+                               VECCOPY(verts[i].x, verts[i].xconst);
+                               VECCOPY(verts[i].v, id->V[i]);
+                       }
+               }
+               else
+               {
+                       VECCOPY(verts[i].txold, id->X[i]);
+                       VECCOPY(verts[i].x, id->X[i]);
+                       VECCOPY(verts[i].v, id->V[i]);
+               }
+       }
+       return 1;
+}
+
+void implicit_set_positions (ClothModifierData *clmd)
+{              
+       Cloth *cloth = clmd->clothObject;
+       ClothVertex *verts = cloth->verts;
+       unsigned int numverts = cloth->numverts, i;
+       Implicit_Data *id = cloth->implicit;
+       
+       for(i = 0; i < numverts; i++)
+       {                               
+               VECCOPY(id->X[i], verts[i].x);
+               VECCOPY(id->V[i], verts[i].v);
+       }
+       if(G.rt > 0)
+               printf("implicit_set_positions\n");     
+}
diff --git a/source/blender/blenkernel/intern/kdop.c b/source/blender/blenkernel/intern/kdop.c
new file mode 100644 (file)
index 0000000..8d2b481
--- /dev/null
@@ -0,0 +1,810 @@
+/*  kdop.c      
+* 
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) Blender Foundation
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): none yet.
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "MEM_guardedalloc.h"
+/* types */
+#include "DNA_curve_types.h"
+#include "DNA_object_types.h"
+#include "DNA_object_force.h"
+#include "DNA_cloth_types.h"   
+#include "DNA_key_types.h"
+#include "DNA_mesh_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_edgehash.h"
+#include "BLI_linklist.h"
+#include "BKE_curve.h"
+#include "BKE_deform.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_cdderivedmesh.h"
+#include "BKE_displist.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+#include "BKE_key.h"
+#include "BKE_mesh.h"
+#include "BKE_object.h"
+#include "BKE_cloth.h"
+#include "BKE_modifier.h"
+#include "BKE_utildefines.h"
+#include "BKE_DerivedMesh.h"
+#include "BIF_editdeform.h"
+#include "BIF_editkey.h"
+#include "DNA_screen_types.h"
+#include "BSE_headerbuttons.h"
+#include "BIF_screen.h"
+#include "BIF_space.h"
+#include "mydevice.h"
+
+
+////////////////////////////////////////////////////////////////////////
+// Additional fastened appending function
+// It uses the link to the last inserted node as start value 
+// for searching the end of the list
+// NEW: in compare to the original function, this one returns
+// the reference to the last inserted node 
+////////////////////////////////////////////////////////////////////////
+LinkNode *BLI_linklist_append_fast(LinkNode **listp, void *ptr) {
+       LinkNode *nlink= MEM_mallocN(sizeof(*nlink), "nlink");
+       LinkNode *node = *listp;
+
+       nlink->link = ptr;
+       nlink->next = NULL;
+
+       if(node == NULL){
+               *listp = nlink;
+       } else {
+               while(node->next != NULL){
+                       node = node->next;   
+               }
+               node->next = nlink;
+       }
+       return nlink;
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+// Bounding Volume Hierarchy Definition
+// 
+// Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
+// Notes: You have to choose the type at compile time ITM
+// Notes: You can choose the tree type --> binary, quad, octree, choose below
+////////////////////////////////////////////////////////////////////////
+
+static float KDOP_AXES[13][3] =
+{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {1, 1, 1}, {1, -1, 1}, {1, 1, -1},
+{1, -1, -1}, {1, 1, 0}, {1, 0, 1}, {0, 1, 1}, {1, -1, 0}, {1, 0, -1},
+{0, 1, -1}
+};
+
+///////////// choose bounding volume here! /////////////
+
+// #define KDOP_26
+
+// #define KDOP_14
+
+// AABB:
+// #define KDOP_8
+
+// OBB: 
+#define KDOP_6
+
+
+
+#ifdef KDOP_26
+#define KDOP_END 13
+#define KDOP_START 0
+#endif
+
+// I didn't test this one!
+#ifdef KDOP_18
+#define KDOP_END 7
+#define KDOP_START 13
+#endif
+
+#ifdef KDOP_14
+#define KDOP_END 7
+#define KDOP_START 0
+#endif
+
+// this is basicly some AABB
+#ifdef KDOP_8
+#define KDOP_END 4
+#define KDOP_START 0
+#endif
+
+// this is basicly some OBB
+#ifdef KDOP_6
+#define KDOP_END 3
+#define KDOP_START 0
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+// Introsort 
+// with permission deriven from the following Java code:
+// http://ralphunden.net/content/tutorials/a-guide-to-introsort/
+// and he derived it from the SUN STL 
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+static int size_threshold = 16;
+/*
+* Common methods for all algorithms
+*/
+DO_INLINE void bvh_exchange(CollisionTree **a, int i, int j)
+{
+       CollisionTree *t=a[i];
+       a[i]=a[j];
+       a[j]=t;
+}
+DO_INLINE int floor_lg(int a)
+{
+       return (int)(floor(log(a)/log(2)));
+}
+
+/*
+* Insertion sort algorithm
+*/
+static void bvh_insertionsort(CollisionTree **a, int lo, int hi, int axis)
+{
+       int i,j;
+       CollisionTree *t;
+       for (i=lo; i < hi; i++)
+       {
+               j=i;
+               t = a[i];
+               while((j!=lo) && (t->bv[axis] < (a[j-1])->bv[axis]))
+               {
+                       a[j] = a[j-1];
+                       j--;
+               }                               
+               a[j] = t;
+       }
+}
+
+static int bvh_partition(CollisionTree **a, int lo, int hi, CollisionTree * x, int axis)
+{
+       int i=lo, j=hi;
+       while (1)
+       {
+               while ((a[i])->bv[axis] < x->bv[axis]) i++;
+               j=j-1;
+               while (x->bv[axis] < (a[j])->bv[axis]) j=j-1;
+               if(!(i < j))
+                       return i;
+               bvh_exchange(a, i,j);
+               i++;
+       }
+}
+
+/*
+* Heapsort algorithm
+*/
+static void bvh_downheap(CollisionTree **a, int i, int n, int lo, int axis)
+{
+       CollisionTree * d = a[lo+i-1];
+       int child;
+       while (i<=n/2)
+       {
+               child = 2*i;
+               if ((child < n) && ((a[lo+child-1])->bv[axis] < (a[lo+child])->bv[axis]))
+               {
+                       child++;
+               }
+               if (!(d->bv[axis] < (a[lo+child-1])->bv[axis])) break;
+               a[lo+i-1] = a[lo+child-1];
+               i = child;
+       }
+       a[lo+i-1] = d;
+}
+
+static void bvh_heapsort(CollisionTree **a, int lo, int hi, int axis)
+{
+       int n = hi-lo, i;
+       for (i=n/2; i>=1; i=i-1)
+       {
+               bvh_downheap(a, i,n,lo, axis);
+       }
+       for (i=n; i>1; i=i-1)
+       {
+               bvh_exchange(a, lo,lo+i-1);
+               bvh_downheap(a, 1,i-1,lo, axis);
+       }
+}
+
+static CollisionTree *bvh_medianof3(CollisionTree **a, int lo, int mid, int hi, int axis) // returns Sortable
+{
+       if ((a[mid])->bv[axis] < (a[lo])->bv[axis])
+       {
+               if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
+                       return a[mid];
+               else
+               {
+                       if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
+                               return a[hi];
+                       else
+                               return a[lo];
+               }
+       }
+       else
+       {
+               if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
+               {
+                       if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
+                               return a[lo];
+                       else
+                               return a[hi];
+               }
+               else
+                       return a[mid];
+       }
+}
+/*
+* Quicksort algorithm modified for Introsort
+*/
+static void bvh_introsort_loop (CollisionTree **a, int lo, int hi, int depth_limit, int axis)
+{
+       int p;
+
+       while (hi-lo > size_threshold)
+       {
+               if (depth_limit == 0)
+               {
+                       bvh_heapsort(a, lo, hi, axis);
+                       return;
+               }
+               depth_limit=depth_limit-1;
+               p=bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo+((hi-lo)/2)+1, hi-1, axis), axis);
+               bvh_introsort_loop(a, p, hi, depth_limit, axis);
+               hi=p;
+       }
+}
+
+DO_INLINE void bvh_sort(CollisionTree **a0, int begin, int end, int axis)
+{
+       if (begin < end)
+       {
+               CollisionTree **a=a0;
+               bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
+               bvh_insertionsort(a, begin, end, axis);
+       }
+}
+DO_INLINE void bvh_sort_along_axis(CollisionTree **face_list, int start, int end, int axis)
+{
+       bvh_sort(face_list, start, end, axis);
+}
+////////////////////////////////////////////////////////////////////////////////////////////////
+void bvh_free(BVH * bvh)
+{
+       LinkNode *search = NULL;
+       CollisionTree *tree = NULL;
+
+       if (bvh) 
+       {
+
+               search = bvh->tree;
+
+               while(search)
+               {
+                       LinkNode *next= search->next;
+                       tree = search->link;
+
+                       MEM_freeN(tree);
+
+                       search = next;
+               }
+
+               BLI_linklist_free(bvh->tree,NULL); 
+               bvh->tree = NULL;
+               
+               if(bvh->current_x)
+                       MEM_freeN(bvh->current_x);
+               if(bvh->current_xold)
+                       MEM_freeN(bvh->current_xold);
+
+               MEM_freeN(bvh);
+               bvh = NULL;
+       }
+}
+
+// only supports x,y,z axis in the moment
+// but we should use a plain and simple function here for speed sake
+DO_INLINE int bvh_largest_axis(float *bv)
+{
+       float middle_point[3];
+
+       middle_point[0] = (bv[1]) - (bv[0]); // x axis
+       middle_point[1] = (bv[3]) - (bv[2]); // y axis
+       middle_point[2] = (bv[5]) - (bv[4]); // z axis
+       if (middle_point[0] > middle_point[1]) 
+       {
+               if (middle_point[0] > middle_point[2])
+                       return 1; // max x axis
+               else
+                       return 5; // max z axis
+       }
+       else 
+       {
+               if (middle_point[1] > middle_point[2])
+                       return 3; // max y axis
+               else
+                       return 5; // max z axis
+       }
+}
+
+// depends on the fact that the BVH's for each face is already build
+DO_INLINE void bvh_calc_DOP_hull_from_faces(BVH * bvh, CollisionTree **tri, int numfaces, float *bv)
+{
+       float newmin,newmax;
+       int i, j;
+       for (j = 0; j < numfaces; j++)
+       {
+               // for all Axes.
+               for (i = KDOP_START; i < KDOP_END; i++)
+               {
+                       newmin = (tri [j])->bv[(2 * i)];        
+                       if ((newmin < bv[(2 * i)]) || (j == 0))
+                       {
+                               bv[(2 * i)] = newmin;
+                       }
+
+                       newmax = (tri [j])->bv[(2 * i) + 1];
+                       if ((newmax > bv[(2 * i) + 1]) || (j == 0))
+                       {
+                               bv[(2 * i) + 1] = newmax;
+                       }
+               }
+       }
+}
+
+DO_INLINE void bvh_calc_DOP_hull_static(BVH * bvh, CollisionTree **tri, int numfaces, float *bv)
+{
+       MFace *tempMFace = bvh->mfaces;
+       float *tempBV = bv;
+       float newminmax;
+       int i, j, k;
+       for (j = 0; j < numfaces; j++)
+       {
+               tempMFace = bvh->mfaces + (tri [j])->tri_index;
+               // 3 or 4 vertices per face.
+               for (k = 0; k < 4; k++)
+               {
+                       int temp = 0;  
+                       // If this is a triangle.
+                       if (k == 3 && !tempMFace->v4)
+                               continue;
+                       // TODO: other name for "temp" this gets all vertices of a face
+                       if (k == 0)
+                               temp = tempMFace->v1;
+                       else if (k == 1)
+                               temp = tempMFace->v2;
+                       else if (k == 2)
+                               temp = tempMFace->v3;
+                       else if (k == 3)
+                               temp = tempMFace->v4;
+                       // for all Axes.
+                       for (i = KDOP_START; i < KDOP_END; i++)
+                       {                               
+                               newminmax = INPR(bvh->current_xold[temp].co, KDOP_AXES[i]);
+                               if ((newminmax < tempBV[(2 * i)]) || (k == 0 && j == 0))
+                                       tempBV[(2 * i)] = newminmax;
+                               if ((newminmax > tempBV[(2 * i) + 1])|| (k == 0 && j == 0))
+                                       tempBV[(2 * i) + 1] = newminmax;
+                       }
+               }
+       }
+}
+
+DO_INLINE void bvh_calc_DOP_hull_moving(BVH * bvh, CollisionTree **tri, int numfaces, float *bv)
+{
+       MFace *tempMFace = bvh->mfaces;
+       float *tempBV = bv;
+       float newminmax;
+       int i, j, k;
+       for (j = 0; j < numfaces; j++)
+       {
+               tempMFace = bvh->mfaces + (tri [j])->tri_index;
+               // 3 or 4 vertices per face.
+               for (k = 0; k < 4; k++)
+               {
+                       int temp = 0;  
+                       // If this is a triangle.
+                       if (k == 3 && !tempMFace->v4)
+                               continue;
+                       // TODO: other name for "temp" this gets all vertices of a face
+                       if (k == 0)
+                               temp = tempMFace->v1;
+                       else if (k == 1)
+                               temp = tempMFace->v2;
+                       else if (k == 2)
+                               temp = tempMFace->v3;
+                       else if (k == 3)
+                               temp = tempMFace->v4;
+                       // for all Axes.
+                       for (i = KDOP_START; i < KDOP_END; i++)
+                       {                               
+                               newminmax = INPR(bvh->current_xold[temp].co, KDOP_AXES[i]);
+                               if ((newminmax < tempBV[(2 * i)]) || (k == 0 && j == 0))
+                                       tempBV[(2 * i)] = newminmax;
+                               if ((newminmax > tempBV[(2 * i) + 1])|| (k == 0 && j == 0))
+                                       tempBV[(2 * i) + 1] = newminmax;
+                               
+                               newminmax = INPR(bvh->current_x[temp].co, KDOP_AXES[i]);
+                               if ((newminmax < tempBV[(2 * i)]) || (k == 0 && j == 0))
+                                       tempBV[(2 * i)] = newminmax;
+                               if ((newminmax > tempBV[(2 * i) + 1])|| (k == 0 && j == 0))
+                                       tempBV[(2 * i) + 1] = newminmax;
+                       }
+               }
+       }
+}
+
+static void bvh_div_env_node(BVH *bvh, CollisionTree *tree, CollisionTree **face_list, unsigned int start, unsigned int end, int lastaxis, LinkNode *nlink)
+{
+       int             i = 0;
+       CollisionTree *newtree = NULL;
+       int laxis = 0, max_nodes=4;
+       unsigned int tstart, tend;
+       LinkNode *nlink1 = nlink;
+       LinkNode *tnlink;
+       tree->traversed = 0;    
+       // Determine which axis to split along
+       laxis = bvh_largest_axis(tree->bv);
+
+       // Sort along longest axis
+       if(laxis!=lastaxis)
+               bvh_sort_along_axis(face_list, start, end, laxis);
+
+       max_nodes = MIN2((end-start + 1 ),4);
+
+       for (i = 0; i < max_nodes; i++)
+       {
+               tree->count_nodes++;
+
+               if(end-start > 4)
+               {
+                       int quarter = ((float)((float)(end - start + 1) / 4.0f));
+                       tstart = start + i * quarter;
+                       tend = tstart + quarter - 1;
+
+                       // be sure that we get all faces
+                       if(i==3)
+                       {
+                               tend = end;
+                       }
+               }
+               else
+               {
+                       tend = tstart = start + i;
+               }
+
+               // Build tree until 4 node left.
+               if ((tend-tstart + 1 ) > 1) 
+               {
+                       newtree = (CollisionTree *)MEM_callocN(sizeof(CollisionTree), "CollisionTree");
+                       tnlink = BLI_linklist_append_fast(&nlink1->next, newtree);
+
+                       newtree->nodes[0] = newtree->nodes[1] = newtree->nodes[2] = newtree->nodes[3] = NULL;
+                       newtree->count_nodes = 0;
+                       newtree->parent = tree;
+                       newtree->isleaf = 0;
+
+                       tree->nodes[i] = newtree;
+
+                       nlink1 = tnlink;
+
+                       bvh_calc_DOP_hull_from_faces(bvh, &face_list[tstart], tend-tstart + 1, tree->nodes[i]->bv);
+
+                       bvh_div_env_node(bvh, tree->nodes[i], face_list, tstart, tend, laxis, nlink1);
+               }
+               else // ok, we have 1 left for this node
+               {
+                       CollisionTree *tnode = face_list[tstart];
+                       tree->nodes[i] = tnode;
+                       tree->nodes[i]->parent = tree;
+               }
+       }
+       return;
+}
+
+/* function cannot be directly called - needs alloced bvh */
+void bvh_build (BVH *bvh)
+{
+       unsigned int i = 0, j = 0, k = 0;
+       CollisionTree **face_list=NULL;
+       CollisionTree *tree=NULL;
+       LinkNode *nlink = NULL;
+       
+       tree = (CollisionTree *)MEM_callocN(sizeof(CollisionTree), "CollisionTree");
+       // TODO: check succesfull alloc
+       BLI_linklist_append(&bvh->tree, tree);
+
+       nlink = bvh->tree;
+
+       if (tree == NULL) 
+       {
+               printf("bvh_build: Out of memory for nodes.\n");
+               bvh_free(bvh);
+               return;
+       }
+       bvh->root = bvh->tree->link;
+       bvh->root->isleaf = 0;
+       bvh->root->parent = NULL;
+       bvh->root->nodes[0] = bvh->root->nodes[1] = bvh->root->nodes[1] = bvh->root->nodes[3] = NULL;
+
+       if(bvh->numfaces<=1)
+       {
+               bvh->root->tri_index = 0;       // Why that? --> only one face there 
+               bvh->root->isleaf = 1;
+               bvh->root->traversed = 0;
+               bvh->root->count_nodes = 0;
+               bvh->leaf_root = bvh->root;
+               bvh->leaf_tree = bvh->root;
+               bvh->root->nextLeaf = NULL;
+               bvh->root->prevLeaf = NULL;
+       }
+       else
+       {
+               // create face boxes            
+               face_list = MEM_callocN (bvh->numfaces * sizeof (CollisionTree *), "CollisionTree");
+               if (face_list == NULL) 
+               {
+                       printf("bvh_build: Out of memory for face_list.\n");
+                       bvh_free(bvh);
+                       return;
+               }
+
+               // create face boxes
+               for(i = 0, k = 0; i < bvh->numfaces; i++)
+               {
+                       LinkNode *tnlink;
+
+                       tree = (CollisionTree *)MEM_callocN(sizeof(CollisionTree), "CollisionTree");
+                       // TODO: check succesfull alloc
+
+                       tnlink = BLI_linklist_append_fast(&nlink->next, tree);
+
+                       face_list[i] = tree;
+                       tree->tri_index = i;
+                       tree->isleaf = 1;
+                       tree->nextLeaf = NULL;
+                       tree->prevLeaf = bvh->leaf_tree;
+                       tree->parent = NULL;
+                       tree->count_nodes = 0;
+
+                       if(i==0)
+                       {
+                               bvh->leaf_tree = bvh->leaf_root = tree;
+                       }
+                       else
+                       {
+                               bvh->leaf_tree->nextLeaf = tree;
+                               bvh->leaf_tree = bvh->leaf_tree->nextLeaf;
+                       }               
+
+                       tree->nodes[0] = tree->nodes[1] = tree->nodes[2] = tree->nodes[3] = NULL;               
+
+                       bvh_calc_DOP_hull_static(bvh, &face_list[i], 1, tree->bv);
+
+                       // inflate the bv with some epsilon
+                       for (j = KDOP_START; j < KDOP_END; j++)
+                       {
+                               tree->bv[(2 * j)] -= bvh->epsilon; // minimum 
+                               tree->bv[(2 * j) + 1] += bvh->epsilon;  // maximum 
+                       }
+                       
+                       nlink = tnlink;
+               }
+               
+               // build root bvh
+               bvh_calc_DOP_hull_from_faces(bvh, face_list, bvh->numfaces, bvh->root->bv);
+
+               // This is the traversal function. 
+               bvh_div_env_node(bvh, bvh->root, face_list, 0, bvh->numfaces-1, 0, nlink);
+               if (face_list)
+                       MEM_freeN(face_list);
+               
+       }
+       
+}
+
+// bvh_overlap - is it possbile for 2 bv's to collide ?
+DO_INLINE int bvh_overlap(float *bv1, float *bv2)
+{
+       int i = 0;
+       for (i = KDOP_START; i < KDOP_END; i++)
+       {
+               // Minimum test.
+               if (bv1[(2 * i)] > bv2[(2 * i) + 1]) 
+               {
+                       return 0;
+               }
+               // Maxiumum test.
+               if (bv2[(2 * i)] > bv1[(2 * i) + 1]) 
+               {
+                       return 0;
+               }
+       }
+       
+       return 1;
+}
+
+/**
+ * bvh_traverse - traverse two bvh trees looking for potential collisions.
+ *
+ * max collisions are n*n collisions --> every triangle collide with
+ * every other triangle that doesn't require any realloc, but uses
+ * much memory
+ */
+int bvh_traverse ( ClothModifierData * clmd, CollisionModifierData * collmd, CollisionTree * tree1, CollisionTree * tree2, float step, CM_COLLISION_RESPONSE collision_response)
+{
+       int i = 0, ret=0;
+       
+       /*
+       // Shouldn't be possible
+       if(!tree1 || !tree2)
+       {
+       printf("Error: no tree there\n");
+       return 0;
+}
+       */      
+       if (bvh_overlap(tree1->bv, tree2->bv)) 
+       {               
+               // Check if this node in the first tree is a leaf
+               if (tree1->isleaf) 
+               {
+                       // Check if this node in the second tree a leaf
+                       if (tree2->isleaf) 
+                       {
+                               // Provide the collision response.
+                               
+                               if(collision_response)
+                                       collision_response (clmd, collmd, tree1, tree2);
+                               return 1;
+                       }
+                       else 
+                       {
+                               // Process the quad tree.
+                               for (i = 0; i < 4; i++)
+                               {
+                                       // Only traverse nodes that exist.
+                                       if (tree2->nodes[i] && bvh_traverse (clmd, collmd, tree1, tree2->nodes[i], step, collision_response))
+                                               ret = 1;
+                               }
+                       }
+               }
+               else 
+               {
+                       // Process the quad tree.
+                       for (i = 0; i < 4; i++)
+                       {
+                               // Only traverse nodes that exist.
+                               if (tree1->nodes [i] && bvh_traverse (clmd, collmd, tree1->nodes[i], tree2, step, collision_response))
+                                       ret = 1;
+                       }
+               }
+       }
+
+       return ret;
+}
+
+// bottom up update of bvh tree:
+// join the 4 children here
+void bvh_join(CollisionTree * tree)
+{
+       int     i = 0, j = 0;
+       if (!tree)
+               return;
+       
+       for (i = 0; i < 4; i++)
+       {
+               if (tree->nodes[i]) 
+               {
+                       for (j = KDOP_START; j < KDOP_END; j++)
+                       {
+                               // update minimum 
+                               if ((tree->nodes[i]->bv[(2 * j)] < tree->bv[(2 * j)]) || (i == 0)) 
+                               {
+                                       tree->bv[(2 * j)] = tree->nodes[i]->bv[(2 * j)];
+                               }
+                               // update maximum 
+                               if ((tree->nodes[i]->bv[(2 * j) + 1] > tree->bv[(2 * j) + 1])|| (i == 0))
+                               {
+                                       tree->bv[(2 * j) + 1] = tree->nodes[i]->bv[(2 * j) + 1];
+                               }
+                       }
+               }
+               else
+                       break;
+       }
+}
+
+// update static bvh
+/* you have to update the bvh position before calling this function */
+void bvh_update(BVH * bvh, int moving)
+{
+       CollisionTree *leaf, *parent;
+       int traversecheck = 1;  // if this is zero we don't go further 
+       unsigned int j = 0;
+       
+       for (leaf = bvh->leaf_root; leaf; leaf = leaf->nextLeaf)
+       {
+               traversecheck = 1;
+               if ((leaf->parent) && (leaf->parent->traversed == leaf->parent->count_nodes)) 
+               {                       
+                       leaf->parent->traversed = 0;
+               }
+               if(!moving)
+                       bvh_calc_DOP_hull_static(bvh, &leaf, 1, leaf->bv);
+               else
+                       bvh_calc_DOP_hull_moving(bvh, &leaf, 1, leaf->bv);
+
+               // inflate the bv with some epsilon
+               for (j = KDOP_START; j < KDOP_END; j++)
+               {                       
+                       leaf->bv[(2 * j)] -= bvh->epsilon; // minimum 
+                       leaf->bv[(2 * j) + 1] += bvh->epsilon;  // maximum 
+               }
+
+               for (parent = leaf->parent; parent; parent = parent->parent)
+               {                       
+                       if (traversecheck) 
+                       {
+                               parent->traversed++;    // we tried to go up in hierarchy 
+                               if (parent->traversed < parent->count_nodes) 
+                               {
+                                       traversecheck = 0;
+
+                                       if (parent->parent) 
+                                       {
+                                               if (parent->parent->traversed == parent->parent->count_nodes) 
+                                               {
+                                                       parent->parent->traversed = 0;
+                                               }
+                                       }
+                                       break;  // we do not need to check further 
+                               }
+                               else 
+                               {
+                                       bvh_join(parent);
+                               }
+                       }
+
+               }
+       }       
+}
+
index 50c6705c6b6f3d7db28b91bc75c1da6685f9f785..d552fdad5c2d3b35a0e5e12d541cf815328c4dd1 100644 (file)
@@ -53,6 +53,7 @@
 #include "MEM_guardedalloc.h"
 
 #include "DNA_armature_types.h"
+#include "DNA_cloth_types.h"
 #include "DNA_effect_types.h"
 #include "DNA_material_types.h"
 #include "DNA_mesh_types.h"
@@ -74,6 +75,7 @@
 #include "BKE_main.h"
 #include "BKE_anim.h"
 #include "BKE_bad_level_calls.h"
+#include "BKE_cloth.h"
 #include "BKE_curve.h"
 #include "BKE_customdata.h"
 #include "BKE_global.h"
@@ -88,6 +90,7 @@
 #include "BKE_object.h"
 #include "BKE_mesh.h"
 #include "BKE_softbody.h"
+#include "BKE_cloth.h"
 #include "BKE_material.h"
 #include "BKE_particle.h"
 #include "BKE_pointcache.h"
@@ -4942,6 +4945,282 @@ static void softbodyModifier_deformVerts(
        sbObjectStep(ob, (float)G.scene->r.cfra, vertexCos, numVerts);
 }
 
+
+/* Cloth */
+
+static void clothModifier_initData(ModifierData *md) 
+{
+       ClothModifierData *clmd = (ClothModifierData*) md;
+       
+       clmd->sim_parms = MEM_callocN(sizeof(SimulationSettings), "cloth sim parms");
+       clmd->coll_parms = MEM_callocN(sizeof(CollisionSettings), "cloth coll parms");
+       
+       /* check for alloc failing */
+       if(!clmd->sim_parms || !clmd->coll_parms)
+               return;
+       
+       cloth_init (clmd);
+       printf("clothModifier_initData\n");
+}
+
+static DerivedMesh *clothModifier_applyModifier(ModifierData *md, Object *ob,
+                DerivedMesh *derivedData, int useRenderParams, int isFinalCalc)
+{
+       ClothModifierData *clmd = (ClothModifierData*) md;
+       DerivedMesh *result=NULL;
+       
+       /* check for alloc failing */
+       if(!clmd->sim_parms || !clmd->coll_parms)
+               return derivedData;
+
+       result = clothModifier_do(clmd, ob, derivedData, useRenderParams, isFinalCalc);
+
+       if(result)
+       {
+               CDDM_calc_normals(result);
+               return result;
+       }
+       return derivedData;
+}
+
+static void clothModifier_updateDepgraph(
+                ModifierData *md, DagForest *forest, Object *ob,
+                DagNode *obNode)
+{
+       ClothModifierData *clmd = (ClothModifierData*) md;
+       
+       Base *base;
+       
+       if(clmd)
+       {
+               for(base = G.scene->base.first; base; base= base->next) 
+               {
+                       Object *ob1= base->object;
+                       if(ob1 != ob)
+                       {
+                               CollisionModifierData *coll_clmd = (CollisionModifierData *)modifiers_findByType(ob1, eModifierType_Collision);
+                               if(coll_clmd)
+                               {
+                                       DagNode *curNode = dag_get_node(forest, ob1);
+                                       dag_add_relation(forest, curNode, obNode, DAG_RL_DATA_DATA|DAG_RL_OB_DATA);
+                               }
+                       }
+               }
+       }       
+}
+
+CustomDataMask clothModifier_requiredDataMask(ModifierData *md)
+{
+       ClothModifierData *clmd = (ClothModifierData *)md;
+       CustomDataMask dataMask = 0;
+
+       /* ask for vertexgroups if we need them */
+       if(clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL)
+               if (clmd->sim_parms->vgroup_mass > 0)
+                       dataMask |= (1 << CD_MDEFORMVERT);
+
+       return dataMask;
+}
+
+static void clothModifier_copyData(ModifierData *md, ModifierData *target)
+{
+       ClothModifierData *clmd = (ClothModifierData*) md;
+       ClothModifierData *tclmd = (ClothModifierData*) target;
+       
+       if(tclmd->sim_parms)
+               MEM_freeN(tclmd->sim_parms);
+       if(tclmd->coll_parms)
+               MEM_freeN(tclmd->coll_parms);   
+       
+       tclmd->sim_parms = MEM_dupallocN(clmd->sim_parms);
+       tclmd->coll_parms = MEM_dupallocN(clmd->coll_parms);
+       
+       tclmd->sim_parms->lastcachedframe = 0;
+}
+
+
+static int clothModifier_dependsOnTime(ModifierData *md)
+{
+       return 1;
+}
+
+static void clothModifier_freeData(ModifierData *md)
+{
+       ClothModifierData *clmd = (ClothModifierData*) md;
+       
+       if (clmd) 
+       {
+               if(G.rt > 0)
+               printf("clothModifier_freeData\n");
+               
+               cloth_free_modifier_extern (clmd);
+               
+               if(clmd->sim_parms)
+                       MEM_freeN(clmd->sim_parms);
+               if(clmd->coll_parms)
+                       MEM_freeN(clmd->coll_parms);    
+       }
+}
+
+/* Collision */
+
+static void collisionModifier_initData(ModifierData *md) 
+{
+       CollisionModifierData *collmd = (CollisionModifierData*) md;
+       
+       collmd->x = NULL;
+       collmd->xnew = NULL;
+       collmd->current_x = NULL;
+       collmd->current_xnew = NULL;
+       collmd->current_v = NULL;
+       collmd->time = -1;
+       collmd->numverts = 0;
+       collmd->tree = NULL;
+}
+
+static void collisionModifier_freeData(ModifierData *md)
+{
+       CollisionModifierData *collmd = (CollisionModifierData*) md;
+       
+       if (collmd) 
+       {
+               if(collmd->tree)
+                       bvh_free(collmd->tree);
+               if(collmd->x)
+                       MEM_freeN(collmd->x);
+               if(collmd->xnew)
+                       MEM_freeN(collmd->xnew);
+               if(collmd->current_x)
+                       MEM_freeN(collmd->current_x);
+               if(collmd->current_xnew)
+                       MEM_freeN(collmd->current_xnew);
+               if(collmd->current_v)
+                       MEM_freeN(collmd->current_v);
+               
+               if(collmd->mfaces)
+                       MEM_freeN(collmd->mfaces);
+               
+               collmd->x = NULL;
+               collmd->xnew = NULL;
+               collmd->current_x = NULL;
+               collmd->current_xnew = NULL;
+               collmd->current_v = NULL;
+               collmd->time = -1;
+               collmd->numverts = 0;
+               collmd->tree = NULL;
+               collmd->mfaces = NULL;
+       }
+}
+
+static int collisionModifier_dependsOnTime(ModifierData *md)
+{
+       return 1;
+}
+
+static void collisionModifier_deformVerts(
+       ModifierData *md, Object *ob, DerivedMesh *derivedData,
+       float (*vertexCos)[3], int numVerts)
+{
+       CollisionModifierData *collmd = (CollisionModifierData*) md;
+       DerivedMesh *dm = NULL;
+       float current_time = 0;
+       unsigned int numverts = 0, i = 0;
+       MVert *tempVert = NULL;
+       
+       /* if possible use/create DerivedMesh */
+       if(derivedData) dm = CDDM_copy(derivedData);
+       else if(ob->type==OB_MESH) dm = CDDM_from_mesh(ob->data, ob);
+       
+       if(!ob->pd)
+       {
+               printf("collisionModifier_deformVerts: Should not happen!\n");
+               return;
+       }
+       
+       if(dm)
+       {
+               CDDM_apply_vert_coords(dm, vertexCos);
+               CDDM_calc_normals(dm);
+               
+               current_time = bsystem_time ( ob, ( float ) G.scene->r.cfra, 0.0 );
+               
+               // printf("current_time %f, collmd->time %f\n", current_time, collmd->time);
+               
+               if(current_time > collmd->time)
+               {       
+                       numverts = dm->getNumVerts ( dm );
+                       
+                       // check if mesh has changed
+                       if(collmd->x && (numverts != collmd->numverts))
+                               collisionModifier_freeData((ModifierData *)collmd);
+                       
+                       if(collmd->time == -1) // first time
+                       {
+                               collmd->x = dm->dupVertArray(dm); // frame start position
+                               
+                               for ( i = 0; i < numverts; i++ )
+                               {
+                                       // we save global positions
+                                       Mat4MulVecfl ( ob->obmat, collmd->x[i].co );
+                               }
+                               
+                               collmd->xnew = MEM_dupallocN(collmd->x); // frame end position
+                               collmd->current_x = MEM_dupallocN(collmd->x); // inter-frame
+                               collmd->current_xnew = MEM_dupallocN(collmd->x); // inter-frame
+                               collmd->current_v = MEM_dupallocN(collmd->x); // inter-frame
+
+                               collmd->numverts = numverts;
+                               
+                               collmd->mfaces = dm->dupFaceArray(dm);
+                               collmd->numfaces = dm->getNumFaces(dm);
+                               
+                               // TODO: epsilon
+                               // create bounding box hierarchy
+                               collmd->tree = bvh_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->x, numverts, ob->pd->pdef_sbift);
+                       }
+                       else if(numverts == collmd->numverts)
+                       {
+                               // put positions to old positions
+                               tempVert = collmd->x;
+                               collmd->x = collmd->xnew;
+                               collmd->xnew = tempVert;
+                               
+                               memcpy(collmd->xnew, dm->getVertArray(dm), numverts*sizeof(MVert));
+                               
+                               for ( i = 0; i < numverts; i++ )
+                               {
+                                       // we save global positions
+                                       Mat4MulVecfl ( ob->obmat, collmd->xnew[i].co );
+                               }
+                               
+                               memcpy(collmd->current_xnew, collmd->x, numverts*sizeof(MVert));
+                               memcpy(collmd->current_x, collmd->x, numverts*sizeof(MVert));
+                               
+                               /* happens on file load (ONLY when i decomment changes in readfile.c */
+                               if(!collmd->tree)
+                               {
+                                       collmd->tree = bvh_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->current_x, numverts, ob->pd->pdef_sbift);
+                               }
+                               else
+                               {
+                                       // recalc static bounding boxes
+                                       bvh_update_from_mvert(collmd->tree, collmd->current_x, numverts, NULL, 0);
+                               }
+                       }
+                       
+                       collmd->time = current_time;
+               }
+               else
+               {
+                       collmd->time = current_time;
+               }
+       }
+       
+       if(dm)
+               dm->release(dm);
+}
+
+
 /* Boolean */
 
 static void booleanModifier_copyData(ModifierData *md, ModifierData *target)
@@ -6765,6 +7044,31 @@ ModifierTypeInfo *modifierType_getInfo(ModifierType type)
                mti->flags = eModifierTypeFlag_AcceptsCVs
                             | eModifierTypeFlag_RequiresOriginalData;
                mti->deformVerts = softbodyModifier_deformVerts;
+       
+               mti = INIT_TYPE(Cloth);
+               mti->type = eModifierTypeType_Nonconstructive;
+               mti->initData = clothModifier_initData;
+               mti->flags = eModifierTypeFlag_AcceptsMesh
+                               | eModifierTypeFlag_RequiresOriginalData;
+                                       // | eModifierTypeFlag_SupportsMapping
+                                       // | eModifierTypeFlag_SupportsEditmode 
+                                       // | eModifierTypeFlag_EnableInEditmode;
+               mti->dependsOnTime = clothModifier_dependsOnTime;
+               mti->freeData = clothModifier_freeData; 
+               mti->requiredDataMask = clothModifier_requiredDataMask;
+               mti->copyData = clothModifier_copyData;
+               mti->applyModifier = clothModifier_applyModifier;
+               mti->updateDepgraph = clothModifier_updateDepgraph;
+               
+               mti = INIT_TYPE(Collision);
+               mti->type = eModifierTypeType_OnlyDeform;
+               mti->initData = collisionModifier_initData;
+               mti->flags = eModifierTypeFlag_AcceptsMesh 
+                               | eModifierTypeFlag_RequiresOriginalData;
+               mti->dependsOnTime = collisionModifier_dependsOnTime;
+               mti->freeData = collisionModifier_freeData; 
+               mti->deformVerts = collisionModifier_deformVerts;
+               // mti->copyData = collisionModifier_copyData;
 
                mti = INIT_TYPE(Boolean);
                mti->type = eModifierTypeType_Nonconstructive;
@@ -7024,6 +7328,13 @@ int modifiers_isSoftbodyEnabled(Object *ob)
        return (md && md->mode & (eModifierMode_Realtime | eModifierMode_Render));
 }
 
+int modifiers_isClothEnabled(Object *ob)
+{
+       ModifierData *md = modifiers_findByType(ob, eModifierType_Cloth);
+
+       return (md && md->mode & (eModifierMode_Realtime | eModifierMode_Render));
+}
+
 int modifiers_isParticleEnabled(Object *ob)
 {
        ModifierData *md = modifiers_findByType(ob, eModifierType_ParticleSystem);
index 3d1a9822ecbc4e7286d7b0082db134e3d73c26d0..f6c5264efc41106e39160daeccb63c3e3a745b1a 100644 (file)
@@ -60,6 +60,7 @@
 #include "DNA_actuator_types.h"
 #include "DNA_brush_types.h"
 #include "DNA_camera_types.h"
+#include "DNA_cloth_types.h"
 #include "DNA_color_types.h"
 #include "DNA_controller_types.h"
 #include "DNA_constraint_types.h"
 
 #include "BKE_action.h"
 #include "BKE_armature.h"
+#include "BKE_cloth.h"
 #include "BKE_colortools.h"
 #include "BKE_constraint.h"
 #include "BKE_curve.h"
@@ -3004,7 +3006,46 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
                        SubsurfModifierData *smd = (SubsurfModifierData*) md;
 
                        smd->emCache = smd->mCache = 0;
-               } else if (md->type==eModifierType_Hook) {
+               }
+               else if (md->type==eModifierType_Cloth) {
+                       ClothModifierData *clmd = (ClothModifierData*) md;
+                       
+                       clmd->clothObject = NULL;
+                       
+                       clmd->sim_parms= newdataadr(fd, clmd->sim_parms);
+                       clmd->coll_parms= newdataadr(fd, clmd->coll_parms);
+                       
+                       clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_LOADED;
+                       clmd->sim_parms->flags &= ~CLOTH_SIMSETTINGS_FLAG_EDITMODE;
+                       
+               }
+               else if (md->type==eModifierType_Collision) {
+                       
+                       CollisionModifierData *collmd = (CollisionModifierData*) md;
+                       /*
+                       // TODO: CollisionModifier should use pointcache 
+                       // + have proper reset events before enabling this
+                       collmd->x = newdataadr(fd, collmd->x);
+                       collmd->xnew = newdataadr(fd, collmd->xnew);
+                       collmd->mfaces = newdataadr(fd, collmd->mfaces);
+                       
+                       collmd->current_x = MEM_callocN(sizeof(MVert)*collmd->numverts,"current_x");
+                       collmd->current_xnew = MEM_callocN(sizeof(MVert)*collmd->numverts,"current_xnew");
+                       collmd->current_v = MEM_callocN(sizeof(MVert)*collmd->numverts,"current_v");
+                       */
+                       
+                       collmd->x = NULL;
+                       collmd->xnew = NULL;
+                       collmd->current_x = NULL;
+                       collmd->current_xnew = NULL;
+                       collmd->current_v = NULL;
+                       collmd->time = -1;
+                       collmd->numverts = 0;
+                       collmd->tree = NULL;
+                       collmd->mfaces = NULL;
+                       
+               }
+               else if (md->type==eModifierType_Hook) {
                        HookModifierData *hmd = (HookModifierData*) md;
 
                        hmd->indexar= newdataadr(fd, hmd->indexar);
@@ -3149,7 +3190,6 @@ static void direct_link_object(FileData *fd, Object *ob)
                sb->bpoint= NULL;       // init pointers so it gets rebuilt nicely
                sb->bspring= NULL;
                sb->scratch= NULL;
-
                /* although not used anymore */
                /* still have to be loaded to be compatible with old files */
                sb->keys= newdataadr(fd, sb->keys);
index f7238e96d8ff2a980ea22211a0dc6ad49d018cc4..8bcc2497dbd0d8c7c17639ac166144f547e026f6 100644 (file)
@@ -107,6 +107,7 @@ Important to know is that 'streaming' has been added to files, for Blender Publi
 #include "DNA_actuator_types.h"
 #include "DNA_brush_types.h"
 #include "DNA_camera_types.h"
+#include "DNA_cloth_types.h"
 #include "DNA_color_types.h"
 #include "DNA_constraint_types.h"
 #include "DNA_controller_types.h"
@@ -155,6 +156,7 @@ Important to know is that 'streaming' has been added to files, for Blender Publi
 #include "BKE_action.h"
 #include "BKE_bad_level_calls.h" // build_seqar (from WHILE_SEQ) free_oops error
 #include "BKE_blender.h"
+#include "BKE_cloth.h"
 #include "BKE_curve.h"
 #include "BKE_customdata.h"
 #include "BKE_constraint.h"
@@ -835,6 +837,24 @@ static void write_modifiers(WriteData *wd, ListBase *modbase)
                        
                        writedata(wd, DATA, sizeof(int)*hmd->totindex, hmd->indexar);
                }
+               else if(md->type==eModifierType_Cloth) {
+                       ClothModifierData *clmd = (ClothModifierData*) md;
+                       
+                       writestruct(wd, DATA, "SimulationSettings", 1, clmd->sim_parms);
+                       writestruct(wd, DATA, "CollisionSettings", 1, clmd->coll_parms);
+                       
+               } 
+               else if (md->type==eModifierType_Collision) {
+                       
+                       CollisionModifierData *collmd = (CollisionModifierData*) md;
+                       /*
+                       // TODO: CollisionModifier should use pointcache 
+                       // + have proper reset events before enabling this
+                       writestruct(wd, DATA, "MVert", collmd->numverts, collmd->x);
+                       writestruct(wd, DATA, "MVert", collmd->numverts, collmd->xnew);
+                       writestruct(wd, DATA, "MFace", collmd->numfaces, collmd->mfaces);
+                       */
+               }
                else if (md->type==eModifierType_MeshDeform) {
                        MeshDeformModifierData *mmd = (MeshDeformModifierData*) md;
                        int size = mmd->dyngridsize;
index a067b3e1489d61e3f501b4ef01995e057c7e1811..28a58f7f9e4e40ac32bffaef615370b103278acf 100644 (file)
@@ -291,6 +291,12 @@ void curvemap_buttons(struct uiBlock *block, struct CurveMapping *cumap, char la
 
 #define B_BAKEABLE_CHANGE              1470
 
+/* Cloth sim button defines */
+#define B_CLOTH_CLEARCACHEALL  1480
+#define B_CLOTH_CLEARCACHEFRAME        1481
+#define B_CLOTH_CHANGEPREROLL  1482
+#define B_CLOTH_RENEW          1483
+
 /* *********************** */
 #define B_WORLDBUTS            1600
 
diff --git a/source/blender/makesdna/DNA_cloth_types.h b/source/blender/makesdna/DNA_cloth_types.h
new file mode 100644 (file)
index 0000000..b7e5478
--- /dev/null
@@ -0,0 +1,655 @@
+/**
+* $Id: DNA_cloth_types.h,v 1.1 2007/08/01 02:28:34 daniel Exp $
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) 2006 by NaN Holding BV.
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel (Genscher)
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+#ifndef DNA_CLOTH_TYPES_H
+#define DNA_CLOTH_TYPES_H
+
+#include "DNA_listBase.h"
+
+/**
+* This struct contains all the global data required to run a simulation.
+* At the time of this writing, this structure contains data appropriate
+* to run a simulation as described in Deformation Constraints in a
+* Mass-Spring Model to Describe Rigid Cloth Behavior by Xavier Provot.
+*
+* I've tried to keep similar, if not exact names for the variables as
+* are presented in the paper.  Where I've changed the concept slightly,
+* as in stepsPerFrame comapred to the time step in the paper, I've used
+* variables with different names to minimize confusion.
+**/
+typedef struct SimulationSettings
+{
+       short   vgroup_mass;    /* optional vertexgroup name for assigning weight.*/
+       short   vgroup_struct;  /* vertex group for scaling structural stiffness */
+       float   mingoal;        /* see SB */
+       int     preroll;        /* How many frames of simulation to do before we start. */
+       float   Cdis;           /* Mechanical damping of springs.               */
+       float   Cvi;            /* Viscous/fluid damping.                       */
+       int     stepsPerFrame;  /* Number of time steps per frame.              */
+       float   gravity [3];    /* Gravity/external force vector.               */
+       float   ufluid [3];     /* Velocity vector of the fluid.                */
+       float   dt;             /* This is the duration of our time step, computed.     */
+       float   mass;           /* The mass of the entire cloth.                */
+       float   structural;     /* Structural spring stiffness.                 */
+       float   shear;          /* Shear spring stiffness.                      */
+       float   bending;        /* Flexion spring stiffness.                    */
+       float   sim_time;
+       int     flags;          /* flags, see CSIMSETT_FLAGS enum above.        */
+       short   solver_type;    /* which solver should be used?         txold   */
+       short   vgroup_bend;    /* vertex group for scaling bending stiffness */
+       float   maxgoal;        /* see SB */
+       float   eff_force_scale;/* Scaling of effector forces (see softbody_calc_forces).*/
+       float   eff_wind_scale; /* Scaling of effector wind (see softbody_calc_forces). */
+       float   sim_time_old;
+       struct  LinkNode *cache; /* UNUSED atm */
+       float   defgoal;
+       int     goalfrict;
+       float   goalspring;
+       int     maxspringlen;   /* in percent!; if tearing enabled, a spring will get cut */
+       int     lastframe;      /* frame on which simulation stops */
+       int     firstframe;     /* frame on which simulation starts */
+       int     lastcachedframe;
+       int     editedframe;    /* which frame is in buffer */
+       int     autoprotect;    /* starting from this frame, cache gets protected  */
+       float   max_bend;       /* max bending scaling value, min is "bending" */
+       float   max_struct;     /* max structural scaling value, min is "structural" */
+       float   max_shear;      /* max shear scaling value, UNUSED */
+       int     firstcachedframe;
+       int pad;
+}
+SimulationSettings;
+
+
+typedef struct CollisionSettings
+{
+       float   epsilon;                /* The radius of a particle in the cloth.               */
+       float   self_friction;          /* Fiction/damping with self contact.                   */
+       float   friction;               /* Friction/damping applied on contact with other object.*/
+       short   collision_type;         /* which collision system is used.                      */
+       short   loop_count;             /* How many iterations for the collision loop.          */
+       struct  LinkNode *collision_list;       /* e.g. pointer to temp memory for collisions */
+       int     flags;                  /* collision flags defined in BKE_cloth.h */
+       float   avg_spring_len;         /* for selfcollision */
+}
+CollisionSettings;
+
+
+/**
+* This structure describes a cloth object against which the
+* simulation can run.
+*
+* The m and n members of this structure represent the assumed
+* rectangular ordered grid for which the original paper is written.
+* At some point they need to disappear and we need to determine out
+* own connectivity of the mesh based on the actual edges in the mesh.
+*
+**/
+typedef struct Cloth
+{
+       struct ClothVertex      *verts;                 /* The vertices that represent this cloth. */
+       struct  LinkNode        *springs;               /* The springs connecting the mesh. */
+       unsigned int            numverts;               /* The number of verts == m * n. */
+       unsigned int            numsprings;             /* The count of springs. */
+       unsigned int            numfaces;
+       unsigned char           old_solver_type;        /* unused, only 1 solver here */
+       unsigned char           pad2;
+       short                   pad3;
+       struct BVH              *tree;                  /* collision tree for this cloth object */
+       struct MFace            *mfaces;
+       struct Implicit_Data    *implicit;              /* our implicit solver connects to this pointer */
+}
+Cloth;
+
+#endif
+/**
+* $Id: DNA_cloth_types.h,v 1.1 2007/08/01 02:28:34 daniel Exp $
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) 2006 by NaN Holding BV.
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel (Genscher)
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+#ifndef DNA_CLOTH_TYPES_H
+#define DNA_CLOTH_TYPES_H
+
+#include "DNA_listBase.h"
+
+/**
+* This struct contains all the global data required to run a simulation.
+* At the time of this writing, this structure contains data appropriate
+* to run a simulation as described in Deformation Constraints in a
+* Mass-Spring Model to Describe Rigid Cloth Behavior by Xavier Provot.
+*
+* I've tried to keep similar, if not exact names for the variables as
+* are presented in the paper.  Where I've changed the concept slightly,
+* as in stepsPerFrame comapred to the time step in the paper, I've used
+* variables with different names to minimize confusion.
+**/
+typedef struct SimulationSettings
+{
+       short   vgroup_mass;    /* optional vertexgroup name for assigning weight.*/
+       short   vgroup_struct;  /* vertex group for scaling structural stiffness */
+       float   mingoal;        /* see SB */
+       int     preroll;        /* How many frames of simulation to do before we start. */
+       float   Cdis;           /* Mechanical damping of springs.               */
+       float   Cvi;            /* Viscous/fluid damping.                       */
+       int     stepsPerFrame;  /* Number of time steps per frame.              */
+       float   gravity [3];    /* Gravity/external force vector.               */
+       float   ufluid [3];     /* Velocity vector of the fluid.                */
+       float   dt;             /* This is the duration of our time step, computed.     */
+       float   mass;           /* The mass of the entire cloth.                */
+       float   structural;     /* Structural spring stiffness.                 */
+       float   shear;          /* Shear spring stiffness.                      */
+       float   bending;        /* Flexion spring stiffness.                    */
+       float   sim_time;
+       int     flags;          /* flags, see CSIMSETT_FLAGS enum above.        */
+       short   solver_type;    /* which solver should be used?         txold   */
+       short   vgroup_bend;    /* vertex group for scaling bending stiffness */
+       float   maxgoal;        /* see SB */
+       float   eff_force_scale;/* Scaling of effector forces (see softbody_calc_forces).*/
+       float   eff_wind_scale; /* Scaling of effector wind (see softbody_calc_forces). */
+       float   sim_time_old;
+       struct  LinkNode *cache; /* UNUSED atm */
+       float   defgoal;
+       int     goalfrict;
+       float   goalspring;
+       int     maxspringlen;   /* in percent!; if tearing enabled, a spring will get cut */
+       int     lastframe;      /* frame on which simulation stops */
+       int     firstframe;     /* frame on which simulation starts */
+       int     lastcachedframe;
+       int     editedframe;    /* which frame is in buffer */
+       int     autoprotect;    /* starting from this frame, cache gets protected  */
+       float   max_bend;       /* max bending scaling value, min is "bending" */
+       float   max_struct;     /* max structural scaling value, min is "structural" */
+       float   max_shear;      /* max shear scaling value, UNUSED */
+       int     firstcachedframe;
+       int pad;
+}
+SimulationSettings;
+
+
+typedef struct CollisionSettings
+{
+       float   epsilon;                /* The radius of a particle in the cloth.               */
+       float   self_friction;          /* Fiction/damping with self contact.                   */
+       float   friction;               /* Friction/damping applied on contact with other object.*/
+       short   collision_type;         /* which collision system is used.                      */
+       short   loop_count;             /* How many iterations for the collision loop.          */
+       struct  LinkNode *collision_list;       /* e.g. pointer to temp memory for collisions */
+       int     flags;                  /* collision flags defined in BKE_cloth.h */
+       float   avg_spring_len;         /* for selfcollision */
+}
+CollisionSettings;
+
+
+/**
+* This structure describes a cloth object against which the
+* simulation can run.
+*
+* The m and n members of this structure represent the assumed
+* rectangular ordered grid for which the original paper is written.
+* At some point they need to disappear and we need to determine out
+* own connectivity of the mesh based on the actual edges in the mesh.
+*
+**/
+typedef struct Cloth
+{
+       struct ClothVertex      *verts;                 /* The vertices that represent this cloth. */
+       struct  LinkNode        *springs;               /* The springs connecting the mesh. */
+       unsigned int            numverts;               /* The number of verts == m * n. */
+       unsigned int            numsprings;             /* The count of springs. */
+       unsigned int            numfaces;
+       unsigned char           old_solver_type;        /* unused, only 1 solver here */
+       unsigned char           pad2;
+       short                   pad3;
+       struct BVH              *tree;                  /* collision tree for this cloth object */
+       struct MFace            *mfaces;
+       struct Implicit_Data    *implicit;              /* our implicit solver connects to this pointer */
+}
+Cloth;
+
+#endif
+/**
+* $Id: DNA_cloth_types.h,v 1.1 2007/08/01 02:28:34 daniel Exp $
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) 2006 by NaN Holding BV.
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel (Genscher)
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+#ifndef DNA_CLOTH_TYPES_H
+#define DNA_CLOTH_TYPES_H
+
+#include "DNA_listBase.h"
+
+/**
+* This struct contains all the global data required to run a simulation.
+* At the time of this writing, this structure contains data appropriate
+* to run a simulation as described in Deformation Constraints in a
+* Mass-Spring Model to Describe Rigid Cloth Behavior by Xavier Provot.
+*
+* I've tried to keep similar, if not exact names for the variables as
+* are presented in the paper.  Where I've changed the concept slightly,
+* as in stepsPerFrame comapred to the time step in the paper, I've used
+* variables with different names to minimize confusion.
+**/
+typedef struct SimulationSettings
+{
+       short   vgroup_mass;    /* optional vertexgroup name for assigning weight.*/
+       short   vgroup_struct;  /* vertex group for scaling structural stiffness */
+       float   mingoal;        /* see SB */
+       int     preroll;        /* How many frames of simulation to do before we start. */
+       float   Cdis;           /* Mechanical damping of springs.               */
+       float   Cvi;            /* Viscous/fluid damping.                       */
+       int     stepsPerFrame;  /* Number of time steps per frame.              */
+       float   gravity [3];    /* Gravity/external force vector.               */
+       float   ufluid [3];     /* Velocity vector of the fluid.                */
+       float   dt;             /* This is the duration of our time step, computed.     */
+       float   mass;           /* The mass of the entire cloth.                */
+       float   structural;     /* Structural spring stiffness.                 */
+       float   shear;          /* Shear spring stiffness.                      */
+       float   bending;        /* Flexion spring stiffness.                    */
+       float   sim_time;
+       int     flags;          /* flags, see CSIMSETT_FLAGS enum above.        */
+       short   solver_type;    /* which solver should be used?         txold   */
+       short   vgroup_bend;    /* vertex group for scaling bending stiffness */
+       float   maxgoal;        /* see SB */
+       float   eff_force_scale;/* Scaling of effector forces (see softbody_calc_forces).*/
+       float   eff_wind_scale; /* Scaling of effector wind (see softbody_calc_forces). */
+       float   sim_time_old;
+       struct  LinkNode *cache; /* UNUSED atm */
+       float   defgoal;
+       int     goalfrict;
+       float   goalspring;
+       int     maxspringlen;   /* in percent!; if tearing enabled, a spring will get cut */
+       int     lastframe;      /* frame on which simulation stops */
+       int     firstframe;     /* frame on which simulation starts */
+       int     lastcachedframe;
+       int     editedframe;    /* which frame is in buffer */
+       int     autoprotect;    /* starting from this frame, cache gets protected  */
+       float   max_bend;       /* max bending scaling value, min is "bending" */
+       float   max_struct;     /* max structural scaling value, min is "structural" */
+       float   max_shear;      /* max shear scaling value, UNUSED */
+       int     firstcachedframe;
+       int pad;
+}
+SimulationSettings;
+
+
+typedef struct CollisionSettings
+{
+       float   epsilon;                /* The radius of a particle in the cloth.               */
+       float   self_friction;          /* Fiction/damping with self contact.                   */
+       float   friction;               /* Friction/damping applied on contact with other object.*/
+       short   collision_type;         /* which collision system is used.                      */
+       short   loop_count;             /* How many iterations for the collision loop.          */
+       struct  LinkNode *collision_list;       /* e.g. pointer to temp memory for collisions */
+       int     flags;                  /* collision flags defined in BKE_cloth.h */
+       float   avg_spring_len;         /* for selfcollision */
+}
+CollisionSettings;
+
+
+/**
+* This structure describes a cloth object against which the
+* simulation can run.
+*
+* The m and n members of this structure represent the assumed
+* rectangular ordered grid for which the original paper is written.
+* At some point they need to disappear and we need to determine out
+* own connectivity of the mesh based on the actual edges in the mesh.
+*
+**/
+typedef struct Cloth
+{
+       struct ClothVertex      *verts;                 /* The vertices that represent this cloth. */
+       struct  LinkNode        *springs;               /* The springs connecting the mesh. */
+       unsigned int            numverts;               /* The number of verts == m * n. */
+       unsigned int            numsprings;             /* The count of springs. */
+       unsigned int            numfaces;
+       unsigned char           old_solver_type;        /* unused, only 1 solver here */
+       unsigned char           pad2;
+       short                   pad3;
+       struct BVH              *tree;                  /* collision tree for this cloth object */
+       struct MFace            *mfaces;
+       struct Implicit_Data    *implicit;              /* our implicit solver connects to this pointer */
+}
+Cloth;
+
+#endif
+/**
+* $Id: DNA_cloth_types.h,v 1.1 2007/08/01 02:28:34 daniel Exp $
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) 2006 by NaN Holding BV.
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel (Genscher)
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+#ifndef DNA_CLOTH_TYPES_H
+#define DNA_CLOTH_TYPES_H
+
+#include "DNA_listBase.h"
+
+/**
+* This struct contains all the global data required to run a simulation.
+* At the time of this writing, this structure contains data appropriate
+* to run a simulation as described in Deformation Constraints in a
+* Mass-Spring Model to Describe Rigid Cloth Behavior by Xavier Provot.
+*
+* I've tried to keep similar, if not exact names for the variables as
+* are presented in the paper.  Where I've changed the concept slightly,
+* as in stepsPerFrame comapred to the time step in the paper, I've used
+* variables with different names to minimize confusion.
+**/
+typedef struct SimulationSettings
+{
+       short   vgroup_mass;    /* optional vertexgroup name for assigning weight.*/
+       short   vgroup_struct;  /* vertex group for scaling structural stiffness */
+       float   mingoal;        /* see SB */
+       int     preroll;        /* How many frames of simulation to do before we start. */
+       float   Cdis;           /* Mechanical damping of springs.               */
+       float   Cvi;            /* Viscous/fluid damping.                       */
+       int     stepsPerFrame;  /* Number of time steps per frame.              */
+       float   gravity [3];    /* Gravity/external force vector.               */
+       float   ufluid [3];     /* Velocity vector of the fluid.                */
+       float   dt;             /* This is the duration of our time step, computed.     */
+       float   mass;           /* The mass of the entire cloth.                */
+       float   structural;     /* Structural spring stiffness.                 */
+       float   shear;          /* Shear spring stiffness.                      */
+       float   bending;        /* Flexion spring stiffness.                    */
+       float   sim_time;
+       int     flags;          /* flags, see CSIMSETT_FLAGS enum above.        */
+       short   solver_type;    /* which solver should be used?         txold   */
+       short   vgroup_bend;    /* vertex group for scaling bending stiffness */
+       float   maxgoal;        /* see SB */
+       float   eff_force_scale;/* Scaling of effector forces (see softbody_calc_forces).*/
+       float   eff_wind_scale; /* Scaling of effector wind (see softbody_calc_forces). */
+       float   sim_time_old;
+       struct  LinkNode *cache; /* UNUSED atm */
+       float   defgoal;
+       int     goalfrict;
+       float   goalspring;
+       int     maxspringlen;   /* in percent!; if tearing enabled, a spring will get cut */
+       int     lastframe;      /* frame on which simulation stops */
+       int     firstframe;     /* frame on which simulation starts */
+       int     lastcachedframe;
+       int     editedframe;    /* which frame is in buffer */
+       int     autoprotect;    /* starting from this frame, cache gets protected  */
+       float   max_bend;       /* max bending scaling value, min is "bending" */
+       float   max_struct;     /* max structural scaling value, min is "structural" */
+       float   max_shear;      /* max shear scaling value, UNUSED */
+       int     firstcachedframe;
+       int pad;
+}
+SimulationSettings;
+
+
+typedef struct CollisionSettings
+{
+       float   epsilon;                /* The radius of a particle in the cloth.               */
+       float   self_friction;          /* Fiction/damping with self contact.                   */
+       float   friction;               /* Friction/damping applied on contact with other object.*/
+       short   collision_type;         /* which collision system is used.                      */
+       short   loop_count;             /* How many iterations for the collision loop.          */
+       struct  LinkNode *collision_list;       /* e.g. pointer to temp memory for collisions */
+       int     flags;                  /* collision flags defined in BKE_cloth.h */
+       float   avg_spring_len;         /* for selfcollision */
+}
+CollisionSettings;
+
+
+/**
+* This structure describes a cloth object against which the
+* simulation can run.
+*
+* The m and n members of this structure represent the assumed
+* rectangular ordered grid for which the original paper is written.
+* At some point they need to disappear and we need to determine out
+* own connectivity of the mesh based on the actual edges in the mesh.
+*
+**/
+typedef struct Cloth
+{
+       struct ClothVertex      *verts;                 /* The vertices that represent this cloth. */
+       struct  LinkNode        *springs;               /* The springs connecting the mesh. */
+       unsigned int            numverts;               /* The number of verts == m * n. */
+       unsigned int            numsprings;             /* The count of springs. */
+       unsigned int            numfaces;
+       unsigned char           old_solver_type;        /* unused, only 1 solver here */
+       unsigned char           pad2;
+       short                   pad3;
+       struct BVH              *tree;                  /* collision tree for this cloth object */
+       struct MFace            *mfaces;
+       struct Implicit_Data    *implicit;              /* our implicit solver connects to this pointer */
+}
+Cloth;
+
+#endif
+/**
+* $Id: DNA_cloth_types.h,v 1.1 2007/08/01 02:28:34 daniel Exp $
+*
+* ***** BEGIN GPL/BL DUAL 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.
+*
+* This program is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program; if not, write to the Free Software Foundation,
+* Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*
+* The Original Code is Copyright (C) 2006 by NaN Holding BV.
+* All rights reserved.
+*
+* The Original Code is: all of this file.
+*
+* Contributor(s): Daniel (Genscher)
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+#ifndef DNA_CLOTH_TYPES_H
+#define DNA_CLOTH_TYPES_H
+
+#include "DNA_listBase.h"
+
+/**
+* This struct contains all the global data required to run a simulation.
+* At the time of this writing, this structure contains data appropriate
+* to run a simulation as described in Deformation Constraints in a
+* Mass-Spring Model to Describe Rigid Cloth Behavior by Xavier Provot.
+*
+* I've tried to keep similar, if not exact names for the variables as
+* are presented in the paper.  Where I've changed the concept slightly,
+* as in stepsPerFrame comapred to the time step in the paper, I've used
+* variables with different names to minimize confusion.
+**/
+typedef struct SimulationSettings
+{
+       short   vgroup_mass;    /* optional vertexgroup name for assigning weight.*/
+       short   vgroup_struct;  /* vertex group for scaling structural stiffness */
+       float   mingoal;        /* see SB */
+       int     preroll;        /* How many frames of simulation to do before we start. */
+       float   Cdis;           /* Mechanical damping of springs.               */
+       float   Cvi;            /* Viscous/fluid damping.                       */
+       int     stepsPerFrame;  /* Number of time steps per frame.              */
+       float   gravity [3];    /* Gravity/external force vector.               */
+       float   ufluid [3];     /* Velocity vector of the fluid.                */
+       float   dt;             /* This is the duration of our time step, computed.     */
+       float   mass;           /* The mass of the entire cloth.                */
+       float   structural;     /* Structural spring stiffness.                 */
+       float   shear;          /* Shear spring stiffness.                      */
+       float   bending;        /* Flexion spring stiffness.                    */
+       float   sim_time;
+       int     flags;          /* flags, see CSIMSETT_FLAGS enum above.        */
+       short   solver_type;    /* which solver should be used?         txold   */
+       short   vgroup_bend;    /* vertex group for scaling bending stiffness */
+       float   maxgoal;        /* see SB */
+       float   eff_force_scale;/* Scaling of effector forces (see softbody_calc_forces).*/
+       float   eff_wind_scale; /* Scaling of effector wind (see softbody_calc_forces). */
+       float   sim_time_old;
+       struct  LinkNode *cache; /* UNUSED atm */
+       float   defgoal;
+       int     goalfrict;
+       float   goalspring;
+       int     maxspringlen;   /* in percent!; if tearing enabled, a spring will get cut */
+       int     lastframe;      /* frame on which simulation stops */
+       int     firstframe;     /* frame on which simulation starts */
+       int     lastcachedframe;
+       int     editedframe;    /* which frame is in buffer */
+       int     autoprotect;    /* starting from this frame, cache gets protected  */
+       float   max_bend;       /* max bending scaling value, min is "bending" */
+       float   max_struct;     /* max structural scaling value, min is "structural" */
+       float   max_shear;      /* max shear scaling value, UNUSED */
+       int     firstcachedframe;
+       int pad;
+}
+SimulationSettings;
+
+
+typedef struct CollisionSettings
+{
+       float   epsilon;                /* The radius of a particle in the cloth.               */
+       float   self_friction;          /* Fiction/damping with self contact.                   */
+       float   friction;               /* Friction/damping applied on contact with other object.*/
+       short   collision_type;         /* which collision system is used.                      */
+       short   loop_count;             /* How many iterations for the collision loop.          */
+       struct  LinkNode *collision_list;       /* e.g. pointer to temp memory for collisions */
+       int     flags;                  /* collision flags defined in BKE_cloth.h */
+       float   avg_spring_len;         /* for selfcollision */
+}
+CollisionSettings;
+
+
+/**
+* This structure describes a cloth object against which the
+* simulation can run.
+*
+* The m and n members of this structure represent the assumed
+* rectangular ordered grid for which the original paper is written.
+* At some point they need to disappear and we need to determine out
+* own connectivity of the mesh based on the actual edges in the mesh.
+*
+**/
+typedef struct Cloth
+{
+       struct ClothVertex      *verts;                 /* The vertices that represent this cloth. */
+       struct  LinkNode        *springs;               /* The springs connecting the mesh. */
+       unsigned int            numverts;               /* The number of verts == m * n. */
+       unsigned int            numsprings;             /* The count of springs. */
+       unsigned int            numfaces;
+       unsigned char           old_solver_type;        /* unused, only 1 solver here */
+       unsigned char           pad2;
+       short                   pad3;
+       struct BVH              *tree;                  /* collision tree for this cloth object */
+       struct MFace            *mfaces;
+       struct Implicit_Data    *implicit;              /* our implicit solver connects to this pointer */
+}
+Cloth;
+
+#endif
index 51f0e9617090afd7b3b93fb28c08d3d6a597705f..1c70508509ba9958f390434f23c9be4d9ca1f414 100644 (file)
@@ -32,6 +32,8 @@ typedef enum ModifierType {
        eModifierType_ParticleSystem,
        eModifierType_ParticleInstance,
        eModifierType_Explode,
+       eModifierType_Cloth,
+       eModifierType_Collision,
        NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -341,6 +343,33 @@ typedef struct SoftbodyModifierData {
        ModifierData modifier;
 } SoftbodyModifierData;
 
+typedef struct ClothModifierData {
+   ModifierData                modifier;
+
+   struct Cloth *clothObject; /* The internal data structure for cloth. */
+   struct SimulationSettings *sim_parms; /* definition is in DNA_cloth_types.h */
+   struct CollisionSettings *coll_parms; /* definition is in DNA_cloth_types.h */
+} ClothModifierData;
+
+typedef struct CollisionModifierData {
+       ModifierData    modifier;
+       
+       struct MVert *x; /* position at the beginning of the frame */
+       struct MVert *xnew; /* position at the end of the frame */
+       struct MVert *xold; /* unsued atm, but was discussed during sprint */
+       struct MVert *current_xnew; /* new position at the actual inter-frame step */
+       struct MVert *current_x; /* position at the actual inter-frame step */
+       struct MVert *current_v; /* position at the actual inter-frame step */
+       
+       struct MFace *mfaces; /* object face data */
+       
+       unsigned int numverts;
+       unsigned int numfaces;
+       int pad;
+       float time;
+       struct BVH *tree;       /* collision tree for this cloth object */
+} CollisionModifierData;
+
 typedef enum {
        eBooleanModifierOp_Intersect,
        eBooleanModifierOp_Union,
index f696c45b3154904c6930b6875c28a90d6ac9056f..ecf2a8e73b0dc483e6cc18ac28aee96878e80940 100644 (file)
@@ -128,6 +128,7 @@ char *includefiles[] = {
        "DNA_brush_types.h",
        "DNA_customdata_types.h",
        "DNA_particle_types.h",
+       "DNA_cloth_types.h",
        // if you add files here, please add them at the end
        // of makesdna.c (this file) as well
 
@@ -1148,4 +1149,5 @@ int main(int argc, char ** argv)
 #include "DNA_brush_types.h"
 #include "DNA_customdata_types.h"
 #include "DNA_particle_types.h"
+#include "DNA_cloth_types.h"
 /* end of list */