Reason of all this work: Commiting my work-in-progress on reviewed collision system (better + general access to kdop, uses less memory, put it into BLI_* namespace and usage defined like existing BLI_kdtree_*). Deleted old kdop.c

24 files changed, 2611 insertions, 1283 deletions

diff --git a/intern/sph/SConscript b/intern/sph/SConscript
new file mode 100644
index 00000000000..52243f767c3
--- /dev/null
+++ b/intern/sph/SConscript
@@ -0,0 +1,11 @@
+#!/usr/bin/python
+import sys
+import os
+Import('env')
+
+sources = env.Glob('intern/*.cpp')
+
+incs = ' . extern intern' 
+defs = ''
+
+env.BlenderLib ('bf_sph', sources, Split(incs), Split(defs), libtype='blender', priority=0 )
diff --git a/intern/sph/extern/sph_extern.h b/intern/sph/extern/sph_extern.h
new file mode 100644
index 00000000000..b4964739212
--- /dev/null
+++ b/intern/sph/extern/sph_extern.h
@@ -0,0 +1,51 @@
+/**
+ * ***** 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): Daniel Genrich.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+
+#ifndef PW_EXTERN_H
+#define PW_EXTERN_H
+
+
+
+#ifdef __cplusplus
+extern "C" { 
+#endif
+/*
+void sph_init_cpp(struct SphModifierData *sphmd);
+void sph_free_cpp(struct SphModifierData *sphmd);
+int sph_simulate_cpp(struct Object *ob, struct SphModifierData *sphmd, float frame, struct ListBase *effectors);
+*/
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif //PW_EXTERN_H	
+
+
diff --git a/intern/sph/intern/sph.cpp b/intern/sph/intern/sph.cpp
new file mode 100644
index 00000000000..f7afa3c34eb
--- /dev/null
+++ b/intern/sph/intern/sph.cpp
@@ -0,0 +1,51 @@
+/*  pw.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.
+*
+* Contributor(s): Daniel Genrich
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include <iostream>
+
+#include "sph_extern.h"
+/*
+extern "C" void sph_init_cpp(struct SphModifierData *sphmd)
+{
+
+}
+/*
+extern "C" void sph_free_cpp(struct SphModifierData *sphmd)
+{
+
+	
+}
+
+extern "C" int sph_simulate_cpp(struct Object *ob, struct SphModifierData *sphmd, float frame, struct ListBase *effectors)
+{
+	
+	return 1;
+}
+*/
diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h
index af920e9762d..f01ed6bbea4 100644
--- a/source/blender/blenkernel/BKE_cloth.h
+++ b/source/blender/blenkernel/BKE_cloth.h
@@ -24,14 +24,14 @@
  *
  * The Original Code is: all of this file.
  *
- * Contributor(s): none yet.
+ * Contributor(s): Daniel Genrich.
  *
  * ***** END GPL LICENSE BLOCK *****
  */
 #ifndef BKE_CLOTH_H
 #define BKE_CLOTH_H
 
-#include "float.h"
+#include <float.h>
 
 #include "BLI_linklist.h"
 #include "BKE_customdata.h"
@@ -49,6 +49,9 @@
 
 #include "BKE_collision.h"
 
+#include "RE_raytrace.h"
+
+
 struct Object;
 
 struct MFace;
@@ -102,7 +105,8 @@ typedef struct Cloth
 	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 BVHTree		*bvhtree;			/* collision tree for this cloth object */
+	struct RayTree 		*selftree;			/* collision tree for this cloth object */
 	struct MFace 		*mfaces;
 	struct Implicit_Data	*implicit; 		/* our implicit solver connects to this pointer */
 	struct Implicit_Data	*implicitEM; 		/* our implicit solver connects to this pointer */
@@ -171,17 +175,10 @@ ClothSpring;
 /* 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_FLAG_AUTOPROTECT = ( 1 << 10 ), /* is autoprotect enabled? */
-	//CLOTH_SIMSETTINGS_FLAG_CCACHE_OUTDATED = (1 << 11),	/* while protected, did cache get outdated? */
 	CLOTH_SIMSETTINGS_FLAG_CCACHE_EDIT = (1 << 12)	/* edit cache in editmode */
 } CLOTH_SIMSETTINGS_FLAGS;
 
@@ -208,6 +205,7 @@ typedef enum
 	CLOTH_SPRING_FLAG_NEEDED = ( 1 << 2 ), // springs has values to be applied
 } CLOTH_SPRINGS_FLAGS;
 
+
 /////////////////////////////////////////////////
 // collision.c
 ////////////////////////////////////////////////
@@ -246,7 +244,7 @@ DerivedMesh *clothModifier_do ( ClothModifierData *clmd,Object *ob, DerivedMesh
 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 );
+void bvhtree_update_from_cloth ( ClothModifierData *clmd, int moving );
 
 // needed for editmesh.c
 void cloth_write_cache ( Object *ob, ClothModifierData *clmd, float framenr );
diff --git a/source/blender/blenkernel/BKE_collision.h b/source/blender/blenkernel/BKE_collision.h
index 7328f9108e3..f0298950f8b 100644
--- a/source/blender/blenkernel/BKE_collision.h
+++ b/source/blender/blenkernel/BKE_collision.h
@@ -32,7 +32,7 @@
 #define BKE_COLLISIONS_H
 
 #include <math.h>
-#include "float.h"
+#include <float.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -47,6 +47,8 @@
 #include "DNA_modifier_types.h"
 #include "DNA_object_types.h"
 
+#include "BLI_kdopbvh.h"
+
 struct Object;
 struct Cloth;
 struct MFace;
@@ -102,10 +104,16 @@ BVH;
 typedef void ( *CM_COLLISION_RESPONSE ) ( ModifierData *md1, ModifierData *md2, CollisionTree *tree1, CollisionTree *tree2 );
 
 // needed for collision.c
-int bvh_traverse ( ModifierData * md1, ModifierData * md2, CollisionTree * tree1, CollisionTree * tree2, float step, CM_COLLISION_RESPONSE collision_response, int selfcollision);
+int bvh_traverse ( ModifierData * md1, ModifierData * md2, CollisionTree * tree1, CollisionTree * tree2, float step, CM_COLLISION_RESPONSE collision_response, int selfcollision );
 
 ////////////////////////////////////////
 
+/* COLLISION FLAGS */
+typedef enum
+{
+	COLLISION_IN_FUTURE = ( 1 << 1 ),
+} COLLISION_FLAGS;
+
 
 ////////////////////////////////////////
 // used for collisions in kdop.c and also collision.c
@@ -119,10 +127,10 @@ typedef struct CollPair
 	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
+	int flag;
 	float time; // collision time, from 0 up to 1
-	unsigned int ap1, ap2, ap3, bp1, bp2, bp3;
-	unsigned int pointsb[4];
+	int ap1, ap2, ap3, bp1, bp2, bp3;
+	int pointsb[4];
 }
 CollPair;
 
@@ -160,8 +168,9 @@ FaceCollPair;
 // 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);
+void bvh_build ( BVH *bvh );
+BVH *bvh_build_from_mvert ( MFace *mfaces, unsigned int numfaces, MVert *x, unsigned int numverts, float epsilon );
+BVHTree *bvhtree_build_from_mvert ( MFace *mfaces, unsigned int numfaces, MVert *x, unsigned int numverts, float epsilon );
 
 // frees the same
 void bvh_free ( BVH * bvh );
@@ -169,20 +178,21 @@ 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->current_xold (static) 
+// update bounding volumes, needs updated positions in  bvh->current_xold (static)
 // and also bvh->current_x if moving==1
-void bvh_update_from_mvert(BVH * bvh, MVert *x, unsigned int numverts, MVert *xnew, int moving);
-void bvh_update(BVH * bvh, int moving);
+void bvh_update_from_mvert ( BVH * bvh, MVert *x, unsigned int numverts, MVert *xnew, int moving );
+void bvh_update ( BVH * bvh, int moving );
+void bvhtree_update_from_mvert ( BVHTree * bvhtree, MFace *faces, int numfaces, MVert *x, MVert *xnew, int numverts, 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);
+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);
+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 );
 
 /////////////////////////////////////////////////
 
diff --git a/source/blender/blenkernel/BKE_effect.h b/source/blender/blenkernel/BKE_effect.h
index 3763a659f2f..15816699285 100644
--- a/source/blender/blenkernel/BKE_effect.h
+++ b/source/blender/blenkernel/BKE_effect.h
@@ -31,6 +31,7 @@
 #ifndef BKE_EFFECT_H
 #define BKE_EFFECT_H
 
+#include "DNA_effect_types.h"
 #include "DNA_object_types.h"
 
 struct Effect;
diff --git a/source/blender/blenkernel/BKE_sph.h b/source/blender/blenkernel/BKE_sph.h
new file mode 100644
index 00000000000..9fa42f5acb2
--- /dev/null
+++ b/source/blender/blenkernel/BKE_sph.h
@@ -0,0 +1,69 @@
+/**
+ * ***** 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): Daniel Genrich.
+ *
+ * ***** END GPL/BL DUAL LICENSE BLOCK *****
+ */
+#ifndef BKE_SPH_H
+#define BKE_SPH_H
+
+
+#include "BKE_DerivedMesh.h"
+#include "BKE_utildefines.h"
+
+#include "BLI_linklist.h"
+
+#include "DNA_modifier_types.h"
+#include "DNA_object_types.h"
+
+void sph_init(SphModifierData *sphmd);
+void sph_free_modifier (SphModifierData *sphmd);
+DerivedMesh *sphModifier_do(SphModifierData *sphmd,Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc);
+int sph_init_all (SphModifierData *sphmd, DerivedMesh *dm, Object *ob);
+
+
+/* SIMULATION FLAGS: goal flags,.. */
+/* These are the bits used in SimSettings.flags. */
+// first 16 (short) flags are used for fluid type identification
+typedef enum
+{
+	SPH_SIMSETTINGS_FLAG_FLUID = ( 1 << 0 ), // Fluid object?
+	SPH_SIMSETTINGS_FLAG_OBSTACLE = ( 1 << 1 ), // Obstacle?
+	SPH_SIMSETTINGS_FLAG_DOMAIN = ( 1 << 2 ), // Fluid domain
+	
+	SPH_SIMSETTINGS_FLAG_GHOSTS = ( 1 << 16 ), // use ghost particles?
+	SPH_SIMSETTINGS_FLAG_OFFLINE = ( 1 << 17 ), // do offline simulation?
+	SPH_SIMSETTINGS_FLAG_MULTIRES = ( 1 << 18 ), // use multires?
+	SPH_SIMSETTINGS_FLAG_VORTICITY = ( 1 << 19 ), // use vorticity enhancement?
+	SPH_SIMSETTINGS_FLAG_BAKING = ( 1 << 20 ), // is domain baking?
+	SPH_SIMSETTINGS_FLAG_INIT = ( 1 << 21 ), // inited?
+} SPH_SIMSETTINGS_FLAGS;
+
+
+#endif //BKE_SPH_H
+
+
+
diff --git a/source/blender/blenkernel/CCGSubSurf.h b/source/blender/blenkernel/CCGSubSurf.h
new file mode 100644
index 00000000000..a8269b7ada0
--- /dev/null
+++ b/source/blender/blenkernel/CCGSubSurf.h
@@ -0,0 +1,152 @@
+/* $Id: CCGSubSurf.h 12931 2007-12-17 18:20:48Z theeth $ */
+
+typedef void* CCGMeshHDL;
+typedef void* CCGVertHDL;
+typedef void* CCGEdgeHDL;
+typedef void* CCGFaceHDL;
+
+typedef struct _CCGVert CCGVert;
+typedef struct _CCGEdge CCGEdge;
+typedef struct _CCGFace CCGFace;
+
+typedef struct _CCGMeshIFC CCGMeshIFC;
+struct _CCGMeshIFC {
+	int			vertUserSize, edgeUserSize, faceUserSize;
+
+	int			vertDataSize;
+};
+
+/***/
+
+typedef void* CCGAllocatorHDL;
+
+typedef struct _CCGAllocatorIFC CCGAllocatorIFC;
+struct _CCGAllocatorIFC {
+	void*		(*alloc)			(CCGAllocatorHDL a, int numBytes);
+	void*		(*realloc)			(CCGAllocatorHDL a, void *ptr, int newSize, int oldSize);
+	void		(*free)				(CCGAllocatorHDL a, void *ptr);
+	void		(*release)			(CCGAllocatorHDL a);
+};
+
+/***/
+
+typedef enum {
+	eCCGError_None = 0,
+
+	eCCGError_InvalidSyncState,
+	eCCGError_InvalidValue,
+} CCGError;
+
+/***/
+
+typedef struct _CCGSubSurf CCGSubSurf;
+
+CCGSubSurf*	ccgSubSurf_new	(CCGMeshIFC *ifc, int subdivisionLevels, CCGAllocatorIFC *allocatorIFC, CCGAllocatorHDL allocator);
+void		ccgSubSurf_free	(CCGSubSurf *ss);
+
+CCGError	ccgSubSurf_sync	(CCGSubSurf *ss);
+
+CCGError	ccgSubSurf_initFullSync		(CCGSubSurf *ss);
+CCGError	ccgSubSurf_initPartialSync	(CCGSubSurf *ss);
+
+CCGError	ccgSubSurf_syncVert		(CCGSubSurf *ss, CCGVertHDL vHDL, void *vertData, int seam, CCGVert **v_r);
+CCGError	ccgSubSurf_syncEdge		(CCGSubSurf *ss, CCGEdgeHDL eHDL, CCGVertHDL e_vHDL0, CCGVertHDL e_vHDL1, float crease, CCGEdge **e_r);
+CCGError	ccgSubSurf_syncFace		(CCGSubSurf *ss, CCGFaceHDL fHDL, int numVerts, CCGVertHDL *vHDLs, CCGFace **f_r);
+
+CCGError	ccgSubSurf_syncVertDel	(CCGSubSurf *ss, CCGVertHDL vHDL);
+CCGError	ccgSubSurf_syncEdgeDel	(CCGSubSurf *ss, CCGEdgeHDL eHDL);
+CCGError	ccgSubSurf_syncFaceDel	(CCGSubSurf *ss, CCGFaceHDL fHDL);
+
+CCGError	ccgSubSurf_processSync	(CCGSubSurf *ss);
+
+CCGError	ccgSubSurf_setSubdivisionLevels		(CCGSubSurf *ss, int subdivisionLevels);
+
+CCGError	ccgSubSurf_setAllowEdgeCreation		(CCGSubSurf *ss, int allowEdgeCreation, float defaultCreaseValue, void *defaultUserData);
+void		ccgSubSurf_getAllowEdgeCreation		(CCGSubSurf *ss, int *allowEdgeCreation_r, float *defaultCreaseValue_r, void *defaultUserData_r);
+
+void		ccgSubSurf_getUseAgeCounts			(CCGSubSurf *ss, int *useAgeCounts_r, int *vertUserOffset_r, int *edgeUserOffset_r, int *faceUserOffset_r);
+CCGError	ccgSubSurf_setUseAgeCounts			(CCGSubSurf *ss, int useAgeCounts, int vertUserOffset, int edgeUserOffset, int faceUserOffset);
+
+CCGError	ccgSubSurf_setCalcVertexNormals		(CCGSubSurf *ss, int useVertNormals, int normalDataOffset);
+
+/***/
+
+int			ccgSubSurf_getNumVerts				(CCGSubSurf *ss);
+int			ccgSubSurf_getNumEdges				(CCGSubSurf *ss);
+int			ccgSubSurf_getNumFaces				(CCGSubSurf *ss);
+
+int			ccgSubSurf_getSubdivisionLevels		(CCGSubSurf *ss);
+int			ccgSubSurf_getEdgeSize				(CCGSubSurf *ss);
+int			ccgSubSurf_getEdgeLevelSize			(CCGSubSurf *ss, int level);
+int			ccgSubSurf_getGridSize				(CCGSubSurf *ss);
+int			ccgSubSurf_getGridLevelSize			(CCGSubSurf *ss, int level);
+
+CCGVert*	ccgSubSurf_getVert					(CCGSubSurf *ss, CCGVertHDL v);
+CCGVertHDL	ccgSubSurf_getVertVertHandle		(CCGSubSurf *ss, CCGVert *v);
+int			ccgSubSurf_getVertNumFaces			(CCGSubSurf *ss, CCGVert *v);
+CCGFace*	ccgSubSurf_getVertFace				(CCGSubSurf *ss, CCGVert *v, int index);
+int			ccgSubSurf_getVertNumEdges			(CCGSubSurf *ss, CCGVert *v);
+CCGEdge*	ccgSubSurf_getVertEdge				(CCGSubSurf *ss, CCGVert *v, int index);
+
+int			ccgSubSurf_getVertAge				(CCGSubSurf *ss, CCGVert *v);
+void*		ccgSubSurf_getVertUserData			(CCGSubSurf *ss, CCGVert *v);
+void*		ccgSubSurf_getVertData				(CCGSubSurf *ss, CCGVert *v);
+void*		ccgSubSurf_getVertLevelData			(CCGSubSurf *ss, CCGVert *v, int level);
+
+CCGEdge*	ccgSubSurf_getEdge					(CCGSubSurf *ss, CCGEdgeHDL e);
+CCGEdgeHDL	ccgSubSurf_getEdgeEdgeHandle		(CCGSubSurf *ss, CCGEdge *e);
+int			ccgSubSurf_getEdgeNumFaces			(CCGSubSurf *ss, CCGEdge *e);
+CCGFace*	ccgSubSurf_getEdgeFace				(CCGSubSurf *ss, CCGEdge *e, int index);
+CCGVert*	ccgSubSurf_getEdgeVert0				(CCGSubSurf *ss, CCGEdge *e);
+CCGVert*	ccgSubSurf_getEdgeVert1				(CCGSubSurf *ss, CCGEdge *e);
+float		ccgSubSurf_getEdgeCrease			(CCGSubSurf *ss, CCGEdge *e);
+
+int			ccgSubSurf_getEdgeAge				(CCGSubSurf *ss, CCGEdge *e);
+void*		ccgSubSurf_getEdgeUserData			(CCGSubSurf *ss, CCGEdge *e);
+void*		ccgSubSurf_getEdgeDataArray			(CCGSubSurf *ss, CCGEdge *e);
+void*		ccgSubSurf_getEdgeData				(CCGSubSurf *ss, CCGEdge *e, int x);
+void*		ccgSubSurf_getEdgeLevelData			(CCGSubSurf *ss, CCGEdge *e, int x, int level);
+
+CCGFace*	ccgSubSurf_getFace					(CCGSubSurf *ss, CCGFaceHDL f);
+CCGFaceHDL	ccgSubSurf_getFaceFaceHandle		(CCGSubSurf *ss, CCGFace *f);
+int			ccgSubSurf_getFaceNumVerts			(CCGSubSurf *ss, CCGFace *f);
+CCGVert*	ccgSubSurf_getFaceVert				(CCGSubSurf *ss, CCGFace *f, int index);
+CCGEdge*	ccgSubSurf_getFaceEdge				(CCGSubSurf *ss, CCGFace *f, int index);
+int			ccgSubSurf_getFaceEdgeIndex			(CCGSubSurf *ss, CCGFace *f, CCGEdge *e);
+
+int			ccgSubSurf_getFaceAge				(CCGSubSurf *ss, CCGFace *f);
+void*		ccgSubSurf_getFaceUserData			(CCGSubSurf *ss, CCGFace *f);
+void*		ccgSubSurf_getFaceCenterData		(CCGSubSurf *ss, CCGFace *f);
+void*		ccgSubSurf_getFaceGridEdgeDataArray	(CCGSubSurf *ss, CCGFace *f, int gridIndex);
+void*		ccgSubSurf_getFaceGridEdgeData		(CCGSubSurf *ss, CCGFace *f, int gridIndex, int x);
+void*		ccgSubSurf_getFaceGridDataArray		(CCGSubSurf *ss, CCGFace *f, int gridIndex);
+void*		ccgSubSurf_getFaceGridData			(CCGSubSurf *ss, CCGFace *f, int gridIndex, int x, int y);
+
+int			ccgSubSurf_getNumFinalVerts		(CCGSubSurf *ss);
+int			ccgSubSurf_getNumFinalEdges		(CCGSubSurf *ss);
+int			ccgSubSurf_getNumFinalFaces		(CCGSubSurf *ss);
+
+/***/
+
+typedef struct _CCGVertIterator CCGVertIterator;
+typedef struct _CCGEdgeIterator CCGEdgeIterator;
+typedef struct _CCGFaceIterator CCGFaceIterator;
+
+CCGVertIterator*	ccgSubSurf_getVertIterator	(CCGSubSurf *ss);
+CCGEdgeIterator*	ccgSubSurf_getEdgeIterator	(CCGSubSurf *ss);
+CCGFaceIterator*	ccgSubSurf_getFaceIterator	(CCGSubSurf *ss);
+
+CCGVert*			ccgVertIterator_getCurrent	(CCGVertIterator *vi);
+int					ccgVertIterator_isStopped	(CCGVertIterator *vi);
+void				ccgVertIterator_next		(CCGVertIterator *vi);
+void				ccgVertIterator_free		(CCGVertIterator *vi);
+
+CCGEdge*			ccgEdgeIterator_getCurrent	(CCGEdgeIterator *ei);
+int					ccgEdgeIterator_isStopped	(CCGEdgeIterator *ei);
+void				ccgEdgeIterator_next		(CCGEdgeIterator *ei);
+void				ccgEdgeIterator_free		(CCGEdgeIterator *ei);
+
+CCGFace*			ccgFaceIterator_getCurrent	(CCGFaceIterator *fi);
+int					ccgFaceIterator_isStopped	(CCGFaceIterator *fi);
+void				ccgFaceIterator_next		(CCGFaceIterator *fi);
+void				ccgFaceIterator_free		(CCGFaceIterator *fi);
diff --git a/source/blender/blenkernel/SConscript b/source/blender/blenkernel/SConscript
index 1bb98239a68..4f77e4f42a7 100644
--- a/source/blender/blenkernel/SConscript
+++ b/source/blender/blenkernel/SConscript
@@ -10,6 +10,7 @@ incs += ' #/intern/iksolver/extern ../blenloader ../quicktime'
 incs += ' #/extern/bullet2/src'
 incs += ' #/intern/bmfont'
 incs += ' #/intern/opennl/extern'
+incs += ' #/intern/sph/extern'
 
 incs += ' ' + env['BF_PYTHON_INC']
 incs += ' ' + env['BF_OPENGL_INC']
diff --git a/source/blender/blenkernel/bmesh_private.h b/source/blender/blenkernel/bmesh_private.h
new file mode 100644
index 00000000000..ad90398bf66
--- /dev/null
+++ b/source/blender/blenkernel/bmesh_private.h
@@ -0,0 +1,71 @@
+/**
+ * BME_private.h    jan 2007
+ *
+ *	low level, 'private' function prototypes for bmesh kernel.
+ *
+ * $Id: BKE_bmesh.h,v 1.00 2007/01/17 17:42:01 Briggs Exp $
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. The Blender
+ * Foundation also sells licenses for use in proprietary software under
+ * the Blender License.  See http://www.blender.org/BL/ for information
+ * about this.	
+ *
+ * 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) 2004 Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): Geoffrey Bantle.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef BMESH_PRIVATE
+#define BMESH_PRIVATE
+
+#include "BKE_bmesh.h"
+
+/*ALLOCATION/DEALLOCATION*/
+struct BME_Vert *BME_addvertlist(struct BME_Mesh *bm, struct BME_Vert *example);
+struct BME_Edge *BME_addedgelist(struct BME_Mesh *bm, struct BME_Vert *v1, struct BME_Vert *v2, struct BME_Edge *example);
+struct BME_Poly *BME_addpolylist(struct BME_Mesh *bm, struct BME_Poly *example); 
+struct BME_Loop *BME_create_loop(struct BME_Mesh *bm, struct BME_Vert *v, struct BME_Edge *e, struct BME_Poly *f, struct BME_Loop *example);
+
+void BME_free_vert(struct BME_Mesh *bm, struct BME_Vert *v);
+void BME_free_edge(struct BME_Mesh *bm, struct BME_Edge *e);
+void BME_free_poly(struct BME_Mesh *bm, struct BME_Poly *f);
+void BME_free_loop(struct BME_Mesh *bm, struct BME_Loop *l);
+void BME_delete_loop(struct BME_Mesh *bm, struct BME_Loop *l);
+
+/*DOUBLE CIRCULAR LINKED LIST FUNCTIONS*/
+void BME_cycle_append(void *h, void *nt);
+int BME_cycle_remove(void *h, void *remn);
+int BME_cycle_validate(int len, void *h);
+/*DISK CYCLE MANAGMENT*/
+int BME_disk_append_edge(struct BME_Edge *e, struct BME_Vert *v);
+void BME_disk_remove_edge(struct BME_Edge *e, struct BME_Vert *v);
+/*RADIAL CYCLE MANAGMENT*/
+void BME_radial_append(struct BME_Edge *e, struct BME_Loop *l);
+void BME_radial_remove_loop(struct BME_Loop *l, struct BME_Edge *e);
+
+/*MISC FUNCTIONS*/
+int BME_edge_swapverts(struct BME_Edge *e, struct BME_Vert *orig, struct BME_Vert *new); /*relink edge*/
+int BME_disk_hasedge(struct BME_Vert *v, struct BME_Edge *e);
+
+/*Error reporting. Shouldnt be called by tools ever.*/
+void BME_error(void);
+#endif
diff --git a/source/blender/blenkernel/intern/cloth.c b/source/blender/blenkernel/intern/cloth.c
index 4be4434dfda..192ebd90faa 100644
--- a/source/blender/blenkernel/intern/cloth.c
+++ b/source/blender/blenkernel/intern/cloth.c
@@ -45,6 +45,8 @@
 
 #include "BKE_pointcache.h"
 
+#include "BLI_kdopbvh.h"
+
 #ifdef _WIN32
 void tstart ( void )
 {}
@@ -151,13 +153,14 @@ void cloth_init ( ClothModifierData *clmd )
 	clmd->sim_parms->goalfrict = 0.0f;
 }
 
-
-BVH *bvh_build_from_cloth (ClothModifierData *clmd, float epsilon)
+BVHTree *bvhtree_build_from_cloth (ClothModifierData *clmd, float epsilon)
 {
-	unsigned int i = 0;
-	BVH	*bvh=NULL;
+	int i;
+	BVHTree *bvhtree;
 	Cloth *cloth = clmd->clothObject;
-	ClothVertex *verts = NULL;
+	ClothVertex *verts;
+	MFace *mfaces;
+	float co[12];
 
 	if(!clmd)
 		return NULL;
@@ -168,69 +171,86 @@ BVH *bvh_build_from_cloth (ClothModifierData *clmd, float epsilon)
 		return NULL;
 	
 	verts = cloth->verts;
+	mfaces = cloth->mfaces;
 	
 	// in the moment, return zero if no faces there
 	if(!cloth->numfaces)
 		return NULL;
 	
-	bvh = MEM_callocN(sizeof(BVH), "BVH");
-	if (bvh == NULL) 
-	{
-		printf("bvh: Out of memory.\n");
-		return NULL;
-	}
-	
-	// springs = cloth->springs;
-	// numsprings = cloth->numsprings;
-
-	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" );
+	// create quadtree with k=26
+	bvhtree = BLI_bvhtree_new(cloth->numfaces, epsilon, 4, 26);
 	
-	if (bvh->current_x == NULL) 
+	// fill tree
+	for(i = 0; i < cloth->numfaces; i++, mfaces++)
 	{
-		printf("bvh: Out of memory.\n");
-		MEM_freeN(bvh);
-		return NULL;
-	}
-	
-	for(i = 0; i < bvh->numverts; i++)
-	{
-		VECCOPY(bvh->current_x[i].co, verts[i].tx);
+		VECCOPY(&co[0*3], verts[mfaces->v1].xold);
+		VECCOPY(&co[1*3], verts[mfaces->v2].xold);
+		VECCOPY(&co[2*3], verts[mfaces->v3].xold);
+		
+		if(mfaces->v4)
+			VECCOPY(&co[3*3], verts[mfaces->v4].xold);
+		
+		BLI_bvhtree_insert(bvhtree, i, co, (mfaces->v4 ? 4 : 3));
 	}
 	
-	bvh_build (bvh);
+	// balance tree
+	BLI_bvhtree_balance(bvhtree);
 	
-	return bvh;
+	return bvhtree;
 }
 
-void bvh_update_from_cloth(ClothModifierData *clmd, int moving)
-{
+void bvhtree_update_from_cloth(ClothModifierData *clmd, int moving)
+{	
 	unsigned int i = 0;
 	Cloth *cloth = clmd->clothObject;
-	BVH *bvh = cloth->tree;
+	BVHTree *bvhtree = cloth->bvhtree;
 	ClothVertex *verts = cloth->verts;
+	MFace *mfaces;
+	float co[12], co_moving[12];
+	int ret = 0;
 	
-	if(!bvh)
+	if(!bvhtree)
 		return;
 	
-	if(cloth->numverts!=bvh->numverts)
-		return;
+	mfaces = cloth->mfaces;
 	
-	if(cloth->verts)
+	// update vertex position in bvh tree
+	if(verts && mfaces)
 	{
-		for(i = 0; i < bvh->numverts; i++)
+		for(i = 0; i < cloth->numfaces; i++, mfaces++)
 		{
-			VECCOPY(bvh->current_x[i].co, verts[i].tx);
-			VECCOPY(bvh->current_xold[i].co, verts[i].txold);
+			VECCOPY(&co[0*3], verts[mfaces->v1].txold);
+			VECCOPY(&co[1*3], verts[mfaces->v2].txold);
+			VECCOPY(&co[2*3], verts[mfaces->v3].txold);
+			
+			if(mfaces->v4)
+				VECCOPY(&co[3*3], verts[mfaces->v4].txold);
+		
+			// copy new locations into array
+			if(moving)
+			{
+				// update moving positions
+				VECCOPY(&co_moving[0*3], verts[mfaces->v1].tx);
+				VECCOPY(&co_moving[1*3], verts[mfaces->v2].tx);
+				VECCOPY(&co_moving[2*3], verts[mfaces->v3].tx);
+				
+				if(mfaces->v4)
+					VECCOPY(&co_moving[3*3], verts[mfaces->v4].tx);
+				
+				ret = BLI_bvhtree_update_node(bvhtree, i, co, co_moving, (mfaces->v4 ? 4 : 3));
+			}
+			else
+			{
+				ret = BLI_bvhtree_update_node(bvhtree, i, co, NULL, (mfaces->v4 ? 4 : 3));
+			}
+			
+			// check if tree is already full
+			if(!ret)
+				break;
 		}
+		
+		BLI_bvhtree_update_tree(bvhtree);
 	}
-	
-	bvh_update(bvh, moving);
 }
 
 int modifiers_indexInObject(Object *ob, ModifierData *md_seek);
@@ -541,8 +561,8 @@ void cloth_free_modifier ( Object *ob, ClothModifierData *clmd )
 		cloth->numsprings = 0;
 
 		// free BVH collision tree
-		if ( cloth->tree )
-			bvh_free ( ( BVH * ) cloth->tree );
+		if ( cloth->bvhtree )
+			BLI_bvhtree_free ( cloth->bvhtree );
 
 		// we save our faces for collision objects
 		if ( cloth->mfaces )
@@ -611,8 +631,8 @@ void cloth_free_modifier_extern ( ClothModifierData *clmd )
 		cloth->numsprings = 0;
 
 		// free BVH collision tree
-		if ( cloth->tree )
-			bvh_free ( ( BVH * ) cloth->tree );
+		if ( cloth->bvhtree )
+			BLI_bvhtree_free ( cloth->bvhtree );
 
 		// we save our faces for collision objects
 		if ( cloth->mfaces )
@@ -810,6 +830,7 @@ static int cloth_from_object(Object *ob, ClothModifierData *clmd, DerivedMesh *d
 		VECCOPY ( verts->xold, verts->x );
 		VECCOPY ( verts->xconst, verts->x );
 		VECCOPY ( verts->txold, verts->x );
+		VECCOPY ( verts->tx, verts->x );
 		VecMulf ( verts->v, 0.0f );
 
 		verts->impulse_count = 0;
@@ -845,12 +866,11 @@ static int cloth_from_object(Object *ob, ClothModifierData *clmd, DerivedMesh *d
 	if(!first)
 		implicit_set_positions(clmd);
 
-	clmd->clothObject->tree = bvh_build_from_cloth ( clmd, clmd->coll_parms->epsilon );
-
+	clmd->clothObject->bvhtree = bvhtree_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 );
diff --git a/source/blender/blenkernel/intern/collision.c b/source/blender/blenkernel/intern/collision.c
index e244ccca306..f3637b4dda2 100644
--- a/source/blender/blenkernel/intern/collision.c
+++ b/source/blender/blenkernel/intern/collision.c
@@ -41,7 +41,6 @@
 #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"
@@ -49,6 +48,38 @@
 
 #include "Bullet-C-Api.h"
 
+#include "BLI_kdopbvh.h"
+#include "BKE_collision.h"
+
+#ifdef _WIN32
+static void start ( void )
+{}
+static void end ( void )
+{
+}
+static double val()
+{
+	return 0;
+}
+#else
+#include <sys/time.h>
+static void mystart ( struct timeval *start, struct timezone *z )
+{
+	gettimeofday ( start, z );
+}
+static void myend ( struct timeval *end, struct timezone *z )
+{
+	gettimeofday ( end,z );
+}
+static double myval ( struct timeval *start, struct timeval *end )
+{
+	double t1, t2;
+	t1 = ( double ) start->tv_sec + ( double ) start->tv_usec/ ( 1000*1000 );
+	t2 = ( double ) end->tv_sec + ( double ) end->tv_usec/ ( 1000*1000 );
+	return t2-t1;
+}
+#endif
+
 /***********************************
 Collision modifier code start
 ***********************************/
@@ -66,58 +97,80 @@ void collision_move_object ( CollisionModifierData *collmd, float step, float pr
 		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 );
 	}
-	bvh_update_from_mvert ( collmd->bvh, collmd->current_x, collmd->numverts, collmd->current_xnew, 1 );
+	bvhtree_update_from_mvert ( collmd->bvhtree, collmd->mfaces, collmd->numfaces, collmd->current_x, collmd->current_xnew, collmd->numverts, 1 );
 }
 
-/* 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 )
+BVHTree *bvhtree_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;
-	}
-
-	// in the moment, return zero if no faces there
-	if ( !numfaces )
-		return NULL;
+	BVHTree *tree;
+	float co[12];
+	int i;
+	MFace *tface = mfaces;
 
-	bvh->epsilon = epsilon;
-	bvh->numfaces = numfaces;
-	bvh->mfaces = mfaces;
+	tree = BLI_bvhtree_new ( numfaces*2, epsilon, 4, 26 );
 
-	// we have no faces, we save seperate points
-	if ( !mfaces )
+	// fill tree
+	for ( i = 0; i < numfaces; i++, tface++ )
 	{
-		bvh->numfaces = numverts;
-	}
+		VECCOPY ( &co[0*3], x[tface->v1].co );
+		VECCOPY ( &co[1*3], x[tface->v2].co );
+		VECCOPY ( &co[2*3], x[tface->v3].co );
+		if ( tface->v4 )
+			VECCOPY ( &co[3*3], x[tface->v4].co );
 
-	bvh->numverts = numverts;
-	bvh->current_x = MEM_dupallocN ( x );
+		BLI_bvhtree_insert ( tree, i, co, ( mfaces->v4 ? 4 : 3 ) );
+	}
 
-	bvh_build ( bvh );
+	// balance tree
+	BLI_bvhtree_balance ( tree );
 
-	return bvh;
+	return tree;
 }
 
-void bvh_update_from_mvert ( BVH * bvh, MVert *x, unsigned int numverts, MVert *xnew, int moving )
+void bvhtree_update_from_mvert ( BVHTree * bvhtree, MFace *faces, int numfaces, MVert *x, MVert *xnew, int numverts, int moving )
 {
-	if ( !bvh )
-		return;
+	int i;
+	MFace *mfaces = faces;
+	float co[12], co_moving[12];
+	int ret = 0;
 
-	if ( numverts!=bvh->numverts )
+	if ( !bvhtree )
 		return;
 
 	if ( x )
-		memcpy ( bvh->current_xold, x, sizeof ( MVert ) * numverts );
+	{
+		for ( i = 0; i < numfaces; i++, mfaces++ )
+		{
+			VECCOPY ( &co[0*3], x[mfaces->v1].co );
+			VECCOPY ( &co[1*3], x[mfaces->v2].co );
+			VECCOPY ( &co[2*3], x[mfaces->v3].co );
+			if ( mfaces->v4 )
+				VECCOPY ( &co[3*3], x[mfaces->v4].co );
+
+			// copy new locations into array
+			if ( moving && xnew )
+			{
+				// update moving positions
+				VECCOPY ( &co_moving[0*3], xnew[mfaces->v1].co );
+				VECCOPY ( &co_moving[1*3], xnew[mfaces->v2].co );
+				VECCOPY ( &co_moving[2*3], xnew[mfaces->v3].co );
+				if ( mfaces->v4 )
+					VECCOPY ( &co_moving[3*3], xnew[mfaces->v4].co );
+
+				ret = BLI_bvhtree_update_node ( bvhtree, i, co, co_moving, ( mfaces->v4 ? 4 : 3 ) );
+			}
+			else
+			{
+				ret = BLI_bvhtree_update_node ( bvhtree, i, co, NULL, ( mfaces->v4 ? 4 : 3 ) );
+			}
 
-	if ( xnew )
-		memcpy ( bvh->current_x, xnew, sizeof ( MVert ) * numverts );
+			// check if tree is already full
+			if ( !ret )
+				break;
+		}
 
-	bvh_update ( bvh, moving );
+		BLI_bvhtree_update_tree ( bvhtree );
+	}
 }
 
 /***********************************
@@ -157,11 +210,11 @@ int gsl_poly_solve_cubic ( float a, float b, float c, float *x0, float *x1, floa
 	else if ( CR2 == CQ3 )
 	{
 		/* this test is actually R2 == Q3, written in a form suitable
-		  for exact computation with integers */
+		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. */
+		considered to be a pair of complex roots z = x +/- epsilon i
+		close to the real axis. */
 
 		float sqrtQ = sqrt ( Q );
 
@@ -419,24 +472,22 @@ DO_INLINE void collision_interpolateOnTriangle ( float to[3], float v1[3], float
 	VECADDMUL ( to, v3, w3 );
 }
 
-int cloth_collision_response_static ( ClothModifierData *clmd, CollisionModifierData *collmd )
+int cloth_collision_response_static ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
 {
 	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;
-	float epsilon2 = collmd->bvh->epsilon;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
 
 	cloth1 = clmd->clothObject;
 
-	search = clmd->coll_parms->collision_list;
-
-	while ( search )
+	for ( ; collpair != collision_end; collpair++ )
 	{
-		collpair = search->link;
+		// only handle static collisions here
+		if ( collpair->flag & COLLISION_IN_FUTURE )
+			continue;
 
 		// compute barycentric coordinates for both collision points
 		collision_compute_barycentric ( collpair->pa,
@@ -530,8 +581,6 @@ int cloth_collision_response_static ( ClothModifierData *clmd, CollisionModifier
 
 			result = 1;
 		}
-
-		search = search->next;
 	}
 
 
@@ -549,51 +598,46 @@ int cloth_collision_response_moving_edges ( ClothModifierData *clmd, ClothModifi
 	return 1;
 }
 
-void cloth_collision_static ( ModifierData *md1, ModifierData *md2, CollisionTree *tree1, CollisionTree *tree2 )
+//Determines collisions on overlap, collisions are writen to collpair[i] and collision+number_collision_found is returned
+CollPair* cloth_collision ( ModifierData *md1, ModifierData *md2, BVHTreeOverlap *overlap, CollPair *collpair )
 {
 	ClothModifierData *clmd = ( ClothModifierData * ) md1;
 	CollisionModifierData *collmd = ( CollisionModifierData * ) md2;
-	CollPair *collpair = NULL;
-	Cloth *cloth1=NULL;
-	MFace *face1=NULL, *face2=NULL;
-	ClothVertex *verts1=NULL;
+	MFace *face1=NULL, *face2 = NULL;
+	ClothVertex *verts1 = clmd->clothObject->verts;
 	double distance = 0;
-	float epsilon = clmd->coll_parms->epsilon;
-	float epsilon2 = ( ( CollisionModifierData * ) md2 )->bvh->epsilon;
-	unsigned int i = 0;
+	float epsilon1 = clmd->coll_parms->epsilon;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
+	int i;
 
+	face1 = & ( clmd->clothObject->mfaces[overlap->indexA] );
+	face2 = & ( collmd->mfaces[overlap->indexB] );
+
+	// check all 4 possible collisions
 	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 )
 		{
+			// fill faceA
 			collpair->ap1 = face1->v1;
 			collpair->ap2 = face1->v2;
 			collpair->ap3 = face1->v3;
 
+			// fill faceB
 			collpair->bp1 = face2->v1;
 			collpair->bp2 = face2->v2;
 			collpair->bp3 = face2->v3;
-
 		}
-
-		if ( i == 1 )
+		else if ( i == 1 )
 		{
 			if ( face1->v4 )
 			{
-				collpair->ap1 = face1->v3;
+				// fill faceA
+				collpair->ap1 = face1->v1;
 				collpair->ap2 = face1->v4;
-				collpair->ap3 = face1->v1;
+				collpair->ap3 = face1->v3;
 
+				// fill faceB
 				collpair->bp1 = face2->v1;
 				collpair->bp2 = face2->v2;
 				collpair->bp3 = face2->v3;
@@ -601,235 +645,215 @@ void cloth_collision_static ( ModifierData *md1, ModifierData *md2, CollisionTre
 			else
 				i++;
 		}
-
 		if ( i == 2 )
 		{
 			if ( face2->v4 )
 			{
+				// fill faceA
 				collpair->ap1 = face1->v1;
 				collpair->ap2 = face1->v2;
 				collpair->ap3 = face1->v3;
 
-				collpair->bp1 = face2->v3;
+				// fill faceB
+				collpair->bp1 = face2->v1;
 				collpair->bp2 = face2->v4;
-				collpair->bp3 = face2->v1;
+				collpair->bp3 = face2->v3;
 			}
 			else
-				i+=2;
+				break;
 		}
-
-		if ( i == 3 )
+		else if ( i == 3 )
 		{
-			if ( ( face1->v4 ) && ( face2->v4 ) )
+			if ( face1->v4 && face2->v4 )
 			{
-				collpair->ap1 = face1->v3;
+				// fill faceA
+				collpair->ap1 = face1->v1;
 				collpair->ap2 = face1->v4;
-				collpair->ap3 = face1->v1;
+				collpair->ap3 = face1->v3;
 
-				collpair->bp1 = face2->v3;
+				// fill faceB
+				collpair->bp1 = face2->v1;
 				collpair->bp2 = face2->v4;
-				collpair->bp3 = face2->v1;
+				collpair->bp3 = face2->v3;
 			}
 			else
-				i++;
+				break;
 		}
 
-		// calc SIPcode (?)
-
-		if ( i < 4 )
-		{
-			// calc distance + normal
 #ifdef WITH_BULLET
-			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 );
+		// 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 );
 #else
-			// just be sure that we don't add anything
-			distance = 2.0 * ( epsilon + epsilon2 + ALMOST_ZERO );
+		// just be sure that we don't add anything
+		distance = 2.0 * ( epsilon1 + epsilon2 + ALMOST_ZERO );
 #endif
-			if ( distance <= ( epsilon + epsilon2 + 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 );
+		if ( distance <= ( epsilon1 + epsilon2 + ALMOST_ZERO ) )
+		{
+			VECCOPY ( collpair->normal, collpair->vector );
+			Normalize ( collpair->normal );
 
-			}
-			else
-			{
-				MEM_freeN ( collpair );
-			}
+			collpair->distance = distance;
+			collpair->flag = 0;
 		}
 		else
 		{
-			MEM_freeN ( collpair );
+			// check for collision in the future
+			collpair->flag |= COLLISION_IN_FUTURE;
 		}
+		collpair++;
 	}
+	return collpair;
 }
 
-int cloth_are_edges_adjacent ( ClothModifierData *clmd, ClothModifierData *coll_clmd, EdgeCollPair *edgecollpair )
+int cloth_are_edges_adjacent ( ClothModifierData *clmd, CollisionModifierData *collmd, EdgeCollPair *edgecollpair )
 {
-	Cloth *cloth1 = NULL, *cloth2 = NULL;
-	ClothVertex *verts1 = NULL, *verts2 = NULL;
+	Cloth *cloth1 = NULL;
+	ClothVertex *verts1 = NULL;
 	float temp[3];
+	MVert *verts2 = collmd->current_x; // old x
 
 	cloth1 = clmd->clothObject;
-	cloth2 = coll_clmd->clothObject;
-
 	verts1 = cloth1->verts;
-	verts2 = cloth2->verts;
 
-	VECSUB ( temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p21].xold );
+	VECSUB ( temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p21].co );
 	if ( ABS ( INPR ( temp, temp ) ) < ALMOST_ZERO )
 		return 1;
 
-	VECSUB ( temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p22].xold );
+	VECSUB ( temp, verts1[edgecollpair->p11].xold, verts2[edgecollpair->p22].co );
 	if ( ABS ( INPR ( temp, temp ) ) < ALMOST_ZERO )
 		return 1;
 
-	VECSUB ( temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p21].xold );
+	VECSUB ( temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p21].co );
 	if ( ABS ( INPR ( temp, temp ) ) < ALMOST_ZERO )
 		return 1;
 
-	VECSUB ( temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p22].xold );
+	VECSUB ( temp, verts1[edgecollpair->p12].xold, verts2[edgecollpair->p22].co );
 	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 )
+void cloth_collision_moving_edges ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair )
 {
 	EdgeCollPair edgecollpair;
-	Cloth *cloth1=NULL, *cloth2=NULL;
-	MFace *face1=NULL, *face2=NULL;
-	ClothVertex *verts1=NULL, *verts2=NULL;
+	Cloth *cloth1=NULL;
+	ClothVertex *verts1=NULL;
 	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];
+	MVert *verts2 = collmd->current_x; // old x
+	MVert *velocity2 = collmd->current_v; // velocity
+	float mintime = 0;
 
 	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++ )
+	for(i = 0; i < 9; i++)
 	{
-		if ( i == 0 )
+		// 9 edge - edge possibilities
+		
+		if(i == 0) // cloth edge: 1-2; coll edge: 1-2
 		{
-			edgecollpair.p11 = face1->v1;
-			edgecollpair.p12 = face1->v2;
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap2;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp2;
 		}
-		else if ( i == 1 )
+		else if(i == 1) // cloth edge: 1-2; coll edge: 2-3
 		{
-			edgecollpair.p11 = face1->v2;
-			edgecollpair.p12 = face1->v3;
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap2;
+			
+			edgecollpair.p21 = collpair->bp2;
+			edgecollpair.p22 = collpair->bp3;
 		}
-		else if ( i == 2 )
+		else if(i == 2) // cloth edge: 1-2; coll edge: 1-3
 		{
-			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
-			}
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap2;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp3;
 		}
-		else if ( i == 3 )
+		else if(i == 3) // cloth edge: 2-3; coll edge: 1-2
 		{
-			if ( face1->v4 )
-			{
-				edgecollpair.p11 = face1->v4;
-				edgecollpair.p12 = face1->v1;
-			}
-			else
-				continue;
+			edgecollpair.p11 = collpair->ap2;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp2;
 		}
-		else
+		else if(i == 4) // cloth edge: 2-3; coll edge: 2-3
 		{
-			edgecollpair.p11 = face1->v3;
-			edgecollpair.p12 = face1->v1;
+			edgecollpair.p11 = collpair->ap2;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp2;
+			edgecollpair.p22 = collpair->bp3;
 		}
-
-
-		for ( j = 0; j < 5; j++ )
+		else if(i == 5) // cloth edge: 2-3; coll edge: 1-3
 		{
-			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 ) )
+			edgecollpair.p11 = collpair->ap2;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp3;
+		}
+		else if(i ==6) // cloth edge: 1-3; coll edge: 1-2
+		{
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp2;
+		}
+		else if(i ==7) // cloth edge: 1-3; coll edge: 2-3
+		{
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp2;
+			edgecollpair.p22 = collpair->bp3;
+		}
+		else if(i == 8) // cloth edge: 1-3; coll edge: 1-3
+		{
+			edgecollpair.p11 = collpair->ap1;
+			edgecollpair.p12 = collpair->ap3;
+			
+			edgecollpair.p21 = collpair->bp1;
+			edgecollpair.p22 = collpair->bp3;
+		}
+		
+		if ( !cloth_are_edges_adjacent ( clmd, collmd, &edgecollpair ) )
+		{
+			// always put coll points in p21/p22
+			VECSUB ( a, verts1[edgecollpair.p12].txold, verts1[edgecollpair.p11].txold );
+			VECSUB ( b, verts1[edgecollpair.p12].tv, verts1[edgecollpair.p11].tv );
+			VECSUB ( c, verts2[edgecollpair.p21].co, verts1[edgecollpair.p11].txold );
+			VECSUB ( d, velocity2[edgecollpair.p21].co, verts1[edgecollpair.p11].tv );
+			VECSUB ( e, verts2[edgecollpair.p22].co, verts1[edgecollpair.p11].txold );
+			VECSUB ( f, velocity2[edgecollpair.p22].co, verts1[edgecollpair.p11].v );
+	
+			numsolutions = cloth_get_collision_time ( a, b, c, d, e, f, solution );
+	
+			for ( k = 0; k < numsolutions; k++ )
 			{
-				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 ) )
 				{
-					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");
-					}
+					//float out_collisionTime = solution[k];
+	
+					// TODO: check for collisions
+	
+					// TODO: put into (edge) collision list
+					
+					mintime = MIN2(mintime, solution[k]);
+	
+					printf("Moving edge found!, mintime: %f\n", mintime);
+					break;
 				}
 			}
 		}
@@ -928,6 +952,7 @@ void cloth_collision_moving_tris ( ClothModifierData *clmd, ClothModifierData *c
 	}
 }
 
+/*
 void cloth_collision_moving ( ClothModifierData *clmd, ClothModifierData *coll_clmd, CollisionTree *tree1, CollisionTree *tree2 )
 {
 	// TODO: check for adjacent
@@ -936,51 +961,68 @@ void cloth_collision_moving ( ClothModifierData *clmd, ClothModifierData *coll_c
 	cloth_collision_moving_tris ( clmd, coll_clmd, tree1, tree2 );
 	cloth_collision_moving_tris ( coll_clmd, clmd, tree2, tree1 );
 }
+*/
 
-void cloth_free_collision_list ( ClothModifierData *clmd )
+int cloth_collision_moving ( ClothModifierData *clmd, CollisionModifierData *collmd, CollPair *collpair, CollPair *collision_end )
 {
-	// free collision list
-	if ( clmd->coll_parms->collision_list )
-	{
-		LinkNode *search = clmd->coll_parms->collision_list;
-		while ( search )
-		{
-			CollPair *coll_pair = search->link;
+	int result = 0;
+	Cloth *cloth1;
+	float w1, w2, w3, u1, u2, u3;
+	float v1[3], v2[3], relativeVelocity[3];
+	float magrelVel;
+	float epsilon2 = BLI_bvhtree_getepsilon ( collmd->bvhtree );
 
-			MEM_freeN ( coll_pair );
-			search = search->next;
-		}
-		BLI_linklist_free ( clmd->coll_parms->collision_list,NULL );
+	cloth1 = clmd->clothObject;
 
-		clmd->coll_parms->collision_list = NULL;
+	for ( ; collpair != collision_end; collpair++ )
+	{
+		// only handle moving collisions here
+		if (!( collpair->flag & COLLISION_IN_FUTURE ))
+			continue;
+		
+		cloth_collision_moving_edges ( clmd, collmd, collpair);
 	}
 }
 
 int cloth_bvh_objcollisions_do ( ClothModifierData * clmd, CollisionModifierData *collmd, float step, float dt )
 {
 	Cloth *cloth = clmd->clothObject;
-	BVH *cloth_bvh= ( BVH * ) cloth->tree;
+	BVHTree *cloth_bvh= ( BVHTree * ) cloth->bvhtree;
 	long i=0, j = 0, numfaces = 0, numverts = 0;
 	ClothVertex *verts = NULL;
+	CollPair *collisions = NULL, *collisions_index = NULL;
 	int ret = 0;
-	unsigned int result = 0;
+	int result = 0;
 	float tnull[3] = {0,0,0};
+	BVHTreeOverlap *overlap = NULL;
+
 
 	numfaces = clmd->clothObject->numfaces;
 	numverts = clmd->clothObject->numverts;
 
 	verts = cloth->verts;
 
-	if ( collmd->bvh )
+	if ( collmd->bvhtree )
 	{
 		/* get pointer to bounding volume hierarchy */
-		BVH *coll_bvh = collmd->bvh;
+		BVHTree *coll_bvh = collmd->bvhtree;
 
 		/* move object to position (step) in time */
 		collision_move_object ( collmd, step + dt, step );
 
 		/* search for overlapping collision pairs */
-		bvh_traverse ( ( ModifierData * ) clmd, ( ModifierData * ) collmd, cloth_bvh->root, coll_bvh->root, step, cloth_collision_static, 0 );
+		overlap = BLI_bvhtree_overlap ( cloth_bvh, coll_bvh, &result );
+
+		collisions = ( CollPair* ) MEM_mallocN ( sizeof ( CollPair ) * result*4, "collision array" ); //*4 since cloth_collision_static can return more than 1 collision
+		collisions_index = collisions;
+
+		for ( i = 0; i < result; i++ )
+		{
+			collisions_index = cloth_collision ( ( ModifierData * ) clmd, ( ModifierData * ) collmd, overlap+i, collisions_index );
+		}
+
+		if ( overlap )
+			MEM_freeN ( overlap );
 	}
 	else
 	{
@@ -994,11 +1036,15 @@ int cloth_bvh_objcollisions_do ( ClothModifierData * clmd, CollisionModifierData
 	{
 		result = 0;
 
-		if ( collmd->bvh )
-			result += cloth_collision_response_static ( clmd, collmd );
+		if ( collmd->bvhtree )
+		{
+			result += cloth_collision_response_static ( clmd, collmd, collisions, collisions_index );
+			result += cloth_collision_moving ( clmd, collmd, collisions, collisions_index );
+		}
 
 		// apply impulses in parallel
 		if ( result )
+		{
 			for ( i = 0; i < numverts; i++ )
 			{
 				// calculate "velocities" (just xnew = xold + v; no dt in v)
@@ -1011,12 +1057,10 @@ int cloth_bvh_objcollisions_do ( ClothModifierData * clmd, CollisionModifierData
 					ret++;
 				}
 			}
-
-		if ( !result )
-			break;
+		}
 	}
 
-	cloth_free_collision_list ( clmd );
+	if ( collisions ) MEM_freeN ( collisions );
 
 	return ret;
 }
@@ -1028,7 +1072,7 @@ int cloth_bvh_objcollision ( ClothModifierData * clmd, float step, float dt )
 	CollisionModifierData *collmd=NULL;
 	Cloth *cloth=NULL;
 	Object *coll_ob=NULL;
-	BVH *cloth_bvh=NULL;
+	BVHTree *cloth_bvh=NULL;
 	long i=0, j = 0, numfaces = 0, numverts = 0;
 	unsigned int result = 0, rounds = 0; // result counts applied collisions; ic is for debug output;
 	ClothVertex *verts = NULL;
@@ -1036,14 +1080,14 @@ int cloth_bvh_objcollision ( ClothModifierData * clmd, float step, float dt )
 	ClothModifierData *tclmd;
 	int collisions = 0, count = 0;
 
-	if ( ( clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ ) || ! ( ( ( Cloth * ) clmd->clothObject )->tree ) )
+	if ( ( clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_COLLOBJ ) || ! ( ( ( Cloth * ) clmd->clothObject )->bvhtree ) )
 	{
 		return 0;
 	}
 
 	cloth = clmd->clothObject;
 	verts = cloth->verts;
-	cloth_bvh = ( BVH * ) cloth->tree;
+	cloth_bvh = ( BVHTree * ) cloth->bvhtree;
 	numfaces = clmd->clothObject->numfaces;
 	numverts = clmd->clothObject->numverts;
 
@@ -1052,12 +1096,11 @@ int cloth_bvh_objcollision ( ClothModifierData * clmd, float step, float dt )
 	////////////////////////////////////////////////////////////
 
 	// update cloth bvh
-	bvh_update_from_cloth ( clmd, 0 ); // 0 means STATIC, 1 means MOVING (see later in this function)
+	bvhtree_update_from_cloth ( clmd, 0 ); // 0 means STATIC, 1 means MOVING (see later in this function)
 
 	do
 	{
 		result = 0;
-		clmd->coll_parms->collision_list = NULL;
 
 		// check all collision objects
 		for ( base = G.scene->base.first; base; base = base->next )
@@ -1126,80 +1169,87 @@ int cloth_bvh_objcollision ( ClothModifierData * clmd, float step, float dt )
 		////////////////////////////////////////////////////////////
 		if ( clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_SELF )
 		{
+
+			MFace *mface = clmd->clothObject->mfaces;
+
 			collisions = 1;
 			verts = cloth->verts; // needed for openMP
 
+
+
+			/*
 			for ( count = 0; count < clmd->coll_parms->self_loop_count; count++ )
 			{
-				if ( collisions )
-				{
-					collisions = 0;
-#pragma omp parallel for private(i,j, collisions) shared(verts, ret)
-					for ( i = 0; i < cloth->numverts; i++ )
-					{
-						for ( j = i + 1; j < cloth->numverts; j++ )
-						{
-							float temp[3];
-							float length = 0;
-							float mindistance = clmd->coll_parms->selfepsilon* ( cloth->verts[i].avg_spring_len + cloth->verts[j].avg_spring_len );
+			if ( collisions )
+			{
+			collisions = 0;
+			#pragma omp parallel for private(i,j, collisions) shared(verts, ret)
+			for ( i = 0; i < cloth->numverts; i++ )
+			{
+			for ( j = i + 1; j < cloth->numverts; j++ )
+			{
+			float temp[3];
+			float length = 0;
+			float mindistance = clmd->coll_parms->selfepsilon* ( cloth->verts[i].avg_spring_len + cloth->verts[j].avg_spring_len );
 
-							if ( clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL )
-							{
-								if ( ( cloth->verts [i].flags & CLOTH_VERT_FLAG_PINNED )
-								        && ( cloth->verts [j].flags & CLOTH_VERT_FLAG_PINNED ) )
-								{
-									continue;
-								}
-							}
+			if ( clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL )
+			{
+			if ( ( cloth->verts [i].flags & CLOTH_VERT_FLAG_PINNED )
+			&& ( cloth->verts [j].flags & CLOTH_VERT_FLAG_PINNED ) )
+			{
+			continue;
+			}
+			}
 
-							VECSUB ( temp, verts[i].tx, verts[j].tx );
+			VECSUB ( temp, verts[i].tx, verts[j].tx );
 
-							if ( ( ABS ( temp[0] ) > mindistance ) || ( ABS ( temp[1] ) > mindistance ) || ( ABS ( temp[2] ) > mindistance ) ) continue;
+			if ( ( ABS ( temp[0] ) > mindistance ) || ( ABS ( temp[1] ) > mindistance ) || ( ABS ( temp[2] ) > mindistance ) ) continue;
 
 							// check for adjacent points (i must be smaller j)
-							if ( BLI_edgehash_haskey ( cloth->edgehash, i, j ) )
-							{
-								continue;
-							}
-
-							length = Normalize ( temp );
-
-							if ( length < mindistance )
-							{
-								float correction = mindistance - length;
-
-								if ( cloth->verts [i].flags & CLOTH_VERT_FLAG_PINNED )
-								{
-									VecMulf ( temp, -correction );
-									VECADD ( verts[j].tx, verts[j].tx, temp );
-								}
-								else if ( cloth->verts [j].flags & CLOTH_VERT_FLAG_PINNED )
-								{
-									VecMulf ( temp, correction );
-									VECADD ( verts[i].tx, verts[i].tx, temp );
-								}
-								else
-								{
-									VecMulf ( temp, -correction*0.5 );
-									VECADD ( verts[j].tx, verts[j].tx, temp );
-
-									VECSUB ( verts[i].tx, verts[i].tx, temp );
-								}
-
-								collisions = 1;
-
-								if ( !ret )
-								{
-#pragma omp critical
-									{
-										ret = 1;
-									}
-								}
-							}
-						}
-					}
-				}
+			if ( BLI_edgehash_haskey ( cloth->edgehash, i, j ) )
+			{
+			continue;
+			}
+
+			length = Normalize ( temp );
+
+			if ( length < mindistance )
+			{
+			float correction = mindistance - length;
+
+			if ( cloth->verts [i].flags & CLOTH_VERT_FLAG_PINNED )
+			{
+			VecMulf ( temp, -correction );
+			VECADD ( verts[j].tx, verts[j].tx, temp );
+			}
+			else if ( cloth->verts [j].flags & CLOTH_VERT_FLAG_PINNED )
+			{
+			VecMulf ( temp, correction );
+			VECADD ( verts[i].tx, verts[i].tx, temp );
+			}
+			else
+			{
+			VecMulf ( temp, -correction*0.5 );
+			VECADD ( verts[j].tx, verts[j].tx, temp );
+
+			VECSUB ( verts[i].tx, verts[i].tx, temp );
+			}
+
+			collisions = 1;
+
+			if ( !ret )
+			{
+			#pragma omp critical
+			{
+			ret = 1;
+			}
+			}
+			}
+			}
+			}
+			}
 			}
+			*/
 			////////////////////////////////////////////////////////////
 
 			////////////////////////////////////////////////////////////
diff --git a/source/blender/blenkernel/intern/kdop.c b/source/blender/blenkernel/intern/kdop.c
deleted file mode 100644
index 3189fe960ad..00000000000
--- a/source/blender/blenkernel/intern/kdop.c
+++ /dev/null
@@ -1,860 +0,0 @@
-/*  kdop.c      
-* 
-*
-* ***** BEGIN GPL LICENSE BLOCK *****
-*
-* This program is free software; you can redistribute it and/or
-* modify it under the terms of the GNU General Public License
-* as published by the Free Software Foundation; either version 2
-* of the License, or (at your option) any later version.
-*
-* 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 LICENSE BLOCK *****
-*/
-
-#include "MEM_guardedalloc.h"
-
-#include "BKE_cloth.h"
-
-#include "DNA_cloth_types.h"
-#include "DNA_mesh_types.h"
-#include "DNA_scene_types.h"
-
-#include "BKE_deform.h"
-#include "BKE_DerivedMesh.h"
-#include "BKE_cdderivedmesh.h"
-#include "BKE_effect.h"
-#include "BKE_global.h"
-#include "BKE_object.h"
-#include "BKE_modifier.h"
-#include "BKE_utildefines.h"
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-
-////////////////////////////////////////////////////////////////////////
-// 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}, {0, 1.0, 0}, {0, 0, 1.0}, {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0},
-{1.0, -1.0, -1.0}, {1.0, 1.0, 0}, {1.0, 0, 1.0}, {0, 1.0, 1.0}, {1.0, -1.0, 0}, {1.0, 0, -1.0},
-{0, 1.0, -1.0}
-};
-
-///////////// choose bounding volume here! /////////////
-
-#define KDOP_26
-
-
-
-#ifdef KDOP_26
-#define KDOP_END 13
-#define KDOP_START 0
-#endif
-
-#ifdef KDOP_18
-#define KDOP_END 13
-#define KDOP_START 7
-#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
-*/
-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
-*/
-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;
-	
-	if(numfaces >0)
-	{
-		// for all Axes.
-		for (i = KDOP_START; i < KDOP_END; i++)
-		{
-			bv[(2 * i)] = (tri [0])->bv[(2 * i)];	
-			bv[(2 * i) + 1] = (tri [0])->bv[(2 * i) + 1];
-		}
-	}
-	
-	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)]))
-			{
-				bv[(2 * i)] = newmin;
-			}
-
-			newmax = (tri [j])->bv[(2 * i) + 1];
-			if ((newmax > bv[(2 * i) + 1]))
-			{
-				bv[(2 * i) + 1] = newmax;
-			}
-		}
-	}
-}
-
-DO_INLINE void bvh_calc_DOP_hull_static(BVH * bvh, CollisionTree **tri, int numfaces, float *bv, CollisionTree *tree)
-{
-	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;
-			}
-		}
-		
-		/* calculate normal of this face */
-		/* (code copied from cdderivedmesh.c) */
-		/*
-		if(tempMFace->v4)
-			CalcNormFloat4(bvh->current_xold[tempMFace->v1].co, bvh->current_xold[tempMFace->v2].co,
-				       bvh->current_xold[tempMFace->v3].co, bvh->current_xold[tempMFace->v4].co, tree->normal);
-		else
-			CalcNormFloat(bvh->current_xold[tempMFace->v1].co, bvh->current_xold[tempMFace->v2].co,
-				      bvh->current_xold[tempMFace->v3].co, tree->normal);
-		
-		tree->alpha = 0;
-		*/
-	}
-}
-
-DO_INLINE void bvh_calc_DOP_hull_moving(BVH * bvh, CollisionTree **tri, int numfaces, float *bv, CollisionTree *tree)
-{
-	MFace *tempMFace = bvh->mfaces;
-	float *tempBV = bv;
-	float newminmax;
-	int i, j, k;
-	
-	/* TODO: calculate normals */
-	
-	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);
-	
-	// maximum is 4 since we have a quad tree
-	max_nodes = MIN2((end-start + 1 ),4);
-
-	for (i = 0; i < max_nodes; i++)
-	{
-		tree->count_nodes++;
-
-		if(end-start+1 > 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;
-	
-	bvh->flags = 0;
-	bvh->leaf_tree = NULL;
-	bvh->leaf_root = NULL;
-	bvh->tree = NULL;
-	
-	if(!bvh->current_x)
-	{
-		bvh_free(bvh);
-		return;
-	}
-	
-	bvh->current_xold = MEM_dupallocN(bvh->current_x);
-	
-	tree = (CollisionTree *)MEM_callocN(sizeof(CollisionTree), "CollisionTree");
-	
-	if (tree == NULL) 
-	{
-		printf("bvh_build: Out of memory for nodes.\n");
-		bvh_free(bvh);
-		return;
-	}
-	
-	BLI_linklist_append(&bvh->tree, tree);
-
-	nlink = bvh->tree;
-
-	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, tree);
-
-			// 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_overlap_self - is it possbile for 2 bv's to selfcollide ?
-DO_INLINE int bvh_overlap_self(CollisionTree * tree1, CollisionTree * tree2)
-{
-	// printf("overlap: %f, q: %f\n", (saacos(INPR(tree1->normal, tree2->normal)) / 2.0)+MAX2(tree1->alpha, tree2->alpha), saacos(INPR(tree1->normal, tree2->normal)));
-	
-	if((saacos(INPR(tree1->normal, tree2->normal)) / 2.0)+MAX2(tree1->alpha, tree2->alpha) > M_PI)
-	{
-		return 1;
-	}
-	else
-		return 0;
-}
-
-/**
- * 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 ( ModifierData * md1, ModifierData * md2, CollisionTree * tree1, CollisionTree * tree2, float step, CM_COLLISION_RESPONSE collision_response, int selfcollision)
-{
-	int i = 0, ret=0, overlap = 0;
-	
-	/*
-	// Shouldn't be possible
-	if(!tree1 || !tree2)
-	{
-	printf("Error: no tree there\n");
-	return 0;
-}
-	*/	
-	
-	if(selfcollision)
-		overlap = bvh_overlap_self(tree1, tree2);
-	else
-		overlap = bvh_overlap(tree1->bv, tree2->bv);
-	
-	if (overlap) 
-	{		
-		// 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 (md1, md2, 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 (md1, md2, tree1, tree2->nodes[i], step, collision_response, selfcollision))
-						ret = 1;
-				}
-			}
-		}
-		else 
-		{
-			// Process the quad tree.
-			for (i = 0; i < 4; i++)
-			{
-				// Only traverse nodes that exist.
-				if (tree1->nodes [i] && bvh_traverse (md1, md2, tree1->nodes[i], tree2, step, collision_response, selfcollision))
-					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;
-	float max = 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];
-				}
-			}
-			
-			/* for selfcollisions */
-			/*
-			if(!i)
-			{
-				tree->alpha = tree->nodes[i]->alpha;
-				VECCOPY(tree->normal, tree->nodes[i]->normal);
-			}
-			else
-			{
-				tree->alpha += saacos(INPR(tree->normal, tree->nodes[i]->normal)) / 2.0;
-				VECADD(tree->normal, tree->normal, tree->nodes[i]->normal);
-				VecMulf(tree->normal, 0.5);
-				max = MAX2(max, tree->nodes[i]->alpha);
-			}
-			*/
-			
-		}
-		else
-			break;
-	}
-	
-	tree->alpha += max;
-}
-
-// 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, leaf);
-		else
-			bvh_calc_DOP_hull_moving(bvh, &leaf, 1, leaf->bv, leaf);
-
-		// 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);
-				}
-			}
-
-		}
-	}	
-}
-
diff --git a/source/blender/blenkernel/intern/modifier.c b/source/blender/blenkernel/intern/modifier.c
index 2a8ba878c41..472fd951168 100644
--- a/source/blender/blenkernel/intern/modifier.c
+++ b/source/blender/blenkernel/intern/modifier.c
@@ -58,6 +58,7 @@
 #include "DNA_material_types.h"
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
+#include "DNA_sph_types.h"
 #include "DNA_modifier_types.h"
 #include "DNA_object_types.h"
 #include "DNA_object_force.h"
@@ -90,6 +91,7 @@
 #include "BKE_object.h"
 #include "BKE_mesh.h"
 #include "BKE_softbody.h"
+#include "BKE_sph.h"
 #include "BKE_cloth.h"
 #include "BKE_material.h"
 #include "BKE_particle.h"
@@ -5178,6 +5180,89 @@ static void clothModifier_freeData(ModifierData *md)
 	}
 }
 
+/* Smooth Particly Hydrodynamics */
+
+static void sphModifier_initData(ModifierData *md) 
+{
+	SphModifierData *sphmd = (SphModifierData*) md;
+	
+	sphmd->sim_parms = MEM_callocN(sizeof(SphSimSettings), "SPH sim parms");
+	sphmd->coll_parms = MEM_callocN(sizeof(SphCollSettings), "SPH coll parms");
+	
+	/* check for alloc failing */
+	if(!sphmd->sim_parms || !sphmd->coll_parms)
+		return;
+	
+	sph_init(sphmd);
+}
+
+static DerivedMesh *sphModifier_applyModifier(ModifierData *md, Object *ob,
+					      DerivedMesh *derivedData, int useRenderParams, int isFinalCalc)
+{
+	SphModifierData *sphmd = (SphModifierData*) md;
+	DerivedMesh *result=NULL;
+	
+	/* check for alloc failing */
+	if(!sphmd->sim_parms || !sphmd->coll_parms)
+		return derivedData;
+
+	result = sphModifier_do(sphmd, ob, derivedData, useRenderParams, isFinalCalc);
+
+	if(result)
+	{
+		CDDM_calc_normals(result);
+		return result;
+	}
+	
+	return derivedData;
+}
+
+static void sphModifier_updateDepgraph(
+				       ModifierData *md, DagForest *forest, Object *ob,
+	   DagNode *obNode)
+{
+	SphModifierData *sphmd = (SphModifierData*) md;
+	
+	Base *base;
+	
+}
+
+CustomDataMask sphModifier_requiredDataMask(ModifierData *md)
+{
+	CustomDataMask dataMask = 0;
+
+	/* ask for vertexgroups if we need them */
+	dataMask |= (1 << CD_MDEFORMVERT);
+
+	return dataMask;
+}
+
+static void sphModifier_copyData(ModifierData *md, ModifierData *target)
+{
+	
+}
+
+
+static int sphModifier_dependsOnTime(ModifierData *md)
+{
+	return 1;
+}
+
+static void sphModifier_freeData(ModifierData *md)
+{
+	SphModifierData *sphmd = (SphModifierData*) md;
+	
+	if (sphmd) 
+	{	
+		sph_free_modifier(sphmd);
+		
+		if(sphmd->sim_parms)
+			MEM_freeN(sphmd->sim_parms);
+		if(sphmd->coll_parms)
+			MEM_freeN(sphmd->coll_parms);
+	}
+}
+
 /* Collision */
 
 static void collisionModifier_initData(ModifierData *md) 
@@ -5191,7 +5276,7 @@ static void collisionModifier_initData(ModifierData *md)
 	collmd->current_v = NULL;
 	collmd->time = -1;
 	collmd->numverts = 0;
-	collmd->bvh = NULL;
+	collmd->bvhtree = NULL;
 }
 
 static void collisionModifier_freeData(ModifierData *md)
@@ -5200,8 +5285,8 @@ static void collisionModifier_freeData(ModifierData *md)
 	
 	if (collmd) 
 	{
-		if(collmd->bvh)
-			bvh_free(collmd->bvh);
+		if(collmd->bvhtree)
+			BLI_bvhtree_free(collmd->bvhtree);
 		if(collmd->x)
 			MEM_freeN(collmd->x);
 		if(collmd->xnew)
@@ -5212,7 +5297,6 @@ static void collisionModifier_freeData(ModifierData *md)
 			MEM_freeN(collmd->current_xnew);
 		if(collmd->current_v)
 			MEM_freeN(collmd->current_v);
-		
 		if(collmd->mfaces)
 			MEM_freeN(collmd->mfaces);
 		
@@ -5223,7 +5307,7 @@ static void collisionModifier_freeData(ModifierData *md)
 		collmd->current_v = NULL;
 		collmd->time = -1;
 		collmd->numverts = 0;
-		collmd->bvh = NULL;
+		collmd->bvhtree = NULL;
 		collmd->mfaces = NULL;
 	}
 }
@@ -5291,9 +5375,8 @@ static void collisionModifier_deformVerts(
 				collmd->mfaces = dm->dupFaceArray(dm);
 				collmd->numfaces = dm->getNumFaces(dm);
 				
-				// TODO: epsilon
 				// create bounding box hierarchy
-				collmd->bvh = bvh_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->x, numverts, ob->pd->pdef_sboft);
+				collmd->bvhtree = bvhtree_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->x, numverts, ob->pd->pdef_sboft);
 				
 				collmd->time = current_time;
 			}
@@ -5316,25 +5399,25 @@ static void collisionModifier_deformVerts(
 				memcpy(collmd->current_x, collmd->x, numverts*sizeof(MVert));
 				
 				/* check if GUI setting has changed for bvh */
-				if(collmd->bvh)
+				if(collmd->bvhtree) 
 				{
-					if(ob->pd->pdef_sboft != collmd->bvh->epsilon)
+					if(ob->pd->pdef_sboft != BLI_bvhtree_getepsilon(collmd->bvhtree))
 					{
-						bvh_free(collmd->bvh);
-						collmd->bvh = bvh_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->current_x, numverts, ob->pd->pdef_sboft);
+						BLI_bvhtree_free(collmd->bvhtree);
+						collmd->bvhtree = bvhtree_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->current_x, numverts, ob->pd->pdef_sboft);
 					}
 			
 				}
 				
-				/* happens on file load (ONLY when i decomment changes in readfile.c */
-				if(!collmd->bvh)
+				/* happens on file load (ONLY when i decomment changes in readfile.c) */
+				if(!collmd->bvhtree)
 				{
-					collmd->bvh = bvh_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->current_x, numverts, ob->pd->pdef_sboft);
+					collmd->bvhtree = bvhtree_build_from_mvert(collmd->mfaces, collmd->numfaces, collmd->current_x, numverts, ob->pd->pdef_sboft);
 				}
 				else
 				{
 					// recalc static bounding boxes
-					bvh_update_from_mvert(collmd->bvh, collmd->current_x, numverts, NULL, 0);
+					bvhtree_update_from_mvert ( collmd->bvhtree, collmd->mfaces, collmd->numfaces, collmd->current_x, NULL, collmd->numverts, 0 );
 				}
 				
 				collmd->time = current_time;
@@ -7106,6 +7189,18 @@ ModifierTypeInfo *modifierType_getInfo(ModifierType type)
 		mti->requiredDataMask = bevelModifier_requiredDataMask;
 		mti->applyModifier = bevelModifier_applyModifier;
 		mti->applyModifierEM = bevelModifier_applyModifierEM;
+		
+		mti = INIT_TYPE(Sph);
+		mti->type = eModifierTypeType_Nonconstructive;
+		mti->flags = eModifierTypeFlag_AcceptsMesh
+				| eModifierTypeFlag_UsesPointCache;
+		mti->initData = sphModifier_initData;
+		mti->copyData = sphModifier_copyData;
+		mti->requiredDataMask = sphModifier_requiredDataMask;
+		mti->applyModifier = sphModifier_applyModifier;
+		mti->dependsOnTime = sphModifier_dependsOnTime;
+		mti->freeData = sphModifier_freeData; 
+		mti->updateDepgraph = sphModifier_updateDepgraph;
 
 		mti = INIT_TYPE(Displace);
 		mti->type = eModifierTypeType_OnlyDeform;
diff --git a/source/blender/blenkernel/intern/sph.c b/source/blender/blenkernel/intern/sph.c
new file mode 100644
index 00000000000..3db1a361d5c
--- /dev/null
+++ b/source/blender/blenkernel/intern/sph.c
@@ -0,0 +1,490 @@
+/*  pw.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.
+*
+* Contributor(s): Daniel Genrich
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include <malloc.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_sph.h"
+
+#include "DNA_sph_types.h"
+
+#include "sph_extern.h"
+
+#include "DNA_effect_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_particle_types.h"
+#include "DNA_scene_types.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+#include "BKE_mesh.h"
+#include "BKE_modifier.h"
+#include "BKE_particle.h"
+
+
+// necessary
+#include "float.h"
+#include "math.h"
+#include "BLI_kdtree.h"
+#include "BLI_arithb.h"
+
+// #include "omp.h"
+
+#ifdef _WIN32
+void ststart ( void )
+{}
+void stend ( void )
+{
+}
+double stval()
+{
+	return 0;
+}
+#else
+#include <sys/time.h>
+
+static struct timeval _ststart, _stend;
+static struct timezone stz;
+
+void ststart(void)
+{
+	gettimeofday(&_ststart, &stz);
+}
+void stend(void)
+{
+	gettimeofday(&_stend,&stz);
+}
+double stval()
+{
+	double t1, t2;
+	t1 =  (double)_ststart.tv_sec + (double)_ststart.tv_usec/(1000*1000);
+	t2 =  (double)_stend.tv_sec + (double)_stend.tv_usec/(1000*1000);
+	return t2-t1;
+}
+#endif 
+
+void sph_init (SphModifierData *sphmd)
+{	
+	/* fill modifier with standard settings */
+	sphmd->sim_parms->timestep = 0.001; // 0.001
+	sphmd->sim_parms->viscosity = 80000.0;
+	sphmd->sim_parms->incompressibility = 8000000.0;
+	sphmd->sim_parms->surfacetension = 8000.0;
+	sphmd->sim_parms->density = 1000.0;
+	sphmd->sim_parms->gravity[2] = -9.81;
+	sphmd->sim_parms->gravity[1] = 0.0;
+	sphmd->sim_parms->gravity[0] = 0.0;
+	sphmd->sim_parms->samplingdistance = 0.02; // length of one cell? 0.01
+	sphmd->sim_parms->smoothinglength = 2.5;
+	sphmd->sim_parms->flags = SPH_SIMSETTINGS_FLAG_GHOSTS | SPH_SIMSETTINGS_FLAG_OFFLINE | SPH_SIMSETTINGS_FLAG_MULTIRES | SPH_SIMSETTINGS_FLAG_DOMAIN;
+	sphmd->sim_parms->computesurfaceevery = 5; // 30000000
+	sphmd->sim_parms->fastmarchingevery = 5;
+	sphmd->sim_parms->dumppovrayevery = 300000;
+	sphmd->sim_parms->dumpimageevery = 30; 
+	sphmd->sim_parms->totaltime = 0.01; // 40.0
+	sphmd->sim_parms->tangentialfriction = 0.1;
+	sphmd->sim_parms->normalfriction = 0.95;
+	sphmd->sim_parms->initiallevel = 1;
+	sphmd->sim_parms->rotation_angle = 0.0;
+	
+	sphmd->sim_parms->rotation_axis[0] = 1.0;
+	sphmd->sim_parms->rotation_axis[1] = 1.0;
+	sphmd->sim_parms->rotation_axis[2] = 1.0;
+	
+	sphmd->sim_parms->rotation_center[0] = 0.0;
+	sphmd->sim_parms->rotation_center[1] = 0.0;
+	sphmd->sim_parms->rotation_center[2] = 0.0;
+
+	sphmd->sim_parms->scenelowerbound[0] = -1.0;
+	sphmd->sim_parms->scenelowerbound[1] = -1.0;
+	sphmd->sim_parms->scenelowerbound[2] = -1.0;
+	
+	sphmd->sim_parms->sceneupperbound[0] = 1.0;
+	sphmd->sim_parms->sceneupperbound[1] = 1.0;
+	sphmd->sim_parms->sceneupperbound[2] = 1.0;
+
+	sphmd->sim_parms->alpha = 2.0;
+	sphmd->sim_parms->beta = 3.0;
+	sphmd->sim_parms->gamma = 1.5;
+	
+	sphmd->sim_parms->numverts = 0;
+	sphmd->sim_parms->numtris = 0;
+	sphmd->sim_parms->verts = NULL;
+	sphmd->sim_parms->tris = NULL;
+	sphmd->sim_parms->normals = NULL;
+	
+	sphmd->sim_parms->resolution = 70;
+	
+	sphmd->sim_parms->co = NULL;
+	sphmd->sim_parms->r = NULL;
+	sphmd->sim_parms->numpart = 0;
+}
+
+void sph_free_modifier (SphModifierData *sphmd)
+{
+	
+	// sph_free_cpp(sphmd);
+	
+	if(sphmd->sim_parms->verts)
+		free(sphmd->sim_parms->verts);
+	if(sphmd->sim_parms->tris)
+		free(sphmd->sim_parms->tris);
+	if(sphmd->sim_parms->normals)
+		free(sphmd->sim_parms->normals);
+	if(sphmd->sim_parms->co)
+		MEM_freeN(sphmd->sim_parms->co);
+	if(sphmd->sim_parms->r)
+		free(sphmd->sim_parms->r);
+	
+}
+
+
+DerivedMesh *sphModifier_do(SphModifierData *sphmd,Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
+{
+	SphSimSettings *sim_parms = sphmd->sim_parms;
+	DerivedMesh *result = NULL;
+	MVert *mvert = NULL;
+	MFace *mface = NULL;
+	int a = 0;
+	float mat[4][4], imat[4][4];
+	Mat4CpyMat4(mat, ob->obmat);
+	Mat4Invert(imat, mat);
+	
+	// only domain is simulated
+	if(!(sim_parms->flags & SPH_SIMSETTINGS_FLAG_DOMAIN) && !(sim_parms->flags & SPH_SIMSETTINGS_FLAG_BAKING))
+	{
+		return dm;
+	}
+	
+	ststart();
+	
+	if(!(sim_parms->flags & SPH_SIMSETTINGS_FLAG_INIT))
+	{
+		if(!sph_init_all (sphmd, dm, ob))
+		{
+			sphmd->sim_parms->flags &= ~SPH_SIMSETTINGS_FLAG_INIT;
+			return dm;
+		}
+	}
+
+	
+	
+	// sph_simulate_cpp(ob, sphmd, 1.0, NULL);
+	
+	stend();
+	
+	printf ( "SPH simulation time: %f\n", ( float ) stval() );
+	
+	if(sim_parms->numverts && sim_parms->numtris)
+	{
+		
+		
+		result = CDDM_new ( sim_parms->numverts, 0, sim_parms->numtris);
+		
+		// copy verts
+		mvert = CDDM_get_verts(result);
+		for(a=0; a<sim_parms->numverts; a++) {
+			MVert *mv = &mvert[a];
+			float *vbCo = &sim_parms->verts[a*3];
+			VECCOPY(mv->co, vbCo);
+			Mat4MulVecfl(imat, mv->co);
+		}
+		
+		mface = CDDM_get_faces(result);
+		for(a=0; a<sim_parms->numtris; a++) {
+			MFace *mf = &mface[a];
+			int *tri = &sim_parms->tris[a*3];
+			mf->v1 = tri[0];
+			mf->v2 = tri[1];
+			mf->v3 = tri[2];
+			test_index_face(mf, NULL, 0, 3);
+		}
+		
+		CDDM_calc_edges ( result );
+		CDDM_calc_normals ( result );
+		
+		return result;
+	}
+	else 
+		return dm;
+}
+
+static void set_min_max(float *min, float *max, float *co)
+{
+	// also calc min + max of bounding box for 3d grid
+	min[0] = MIN2(min[0], co[0]);
+	min[1] = MIN2(min[1], co[1]);
+	min[2] = MIN2(min[2], co[2]);
+		
+	max[0] = MAX2(max[0], co[0]);
+	max[1] = MAX2(max[1], co[1]);
+	max[2] = MAX2(max[2], co[2]);
+}
+
+long calc_distance_field(Object *ob, DerivedMesh *dm, SphModifierData *sphmd, float mat[4][4])
+{
+	int numverts = dm->getNumVerts(dm);
+	int numfaces = dm->getNumFaces(dm);
+	MVert *mvert = dm->getVertArray(dm);
+	MFace *mface = dm->getFaceArray(dm);
+	int i, j, k;
+	KDTree *tree;
+	float co[3], co1[3], co2[3], co3[3], co4[3];
+	int index;
+	float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX}, slice, maxdir;
+	int resx, resy, resz;
+	int maxres = 20;
+	float *dist;
+	float *normals;
+	int totpart = 0;
+	float *cos = NULL;
+	int maxpart = 0;
+	
+	printf("calc_distance_field\n");
+	
+	dist = MEM_callocN(maxres*maxres*maxres*sizeof(float), "distance_field");
+	normals = MEM_callocN(numfaces*3*sizeof(float), "triangle_normals");
+	
+	/////////////////////////////////////////////////
+	// create + fill + balance kdtree
+	/////////////////////////////////////////////////
+	tree = BLI_kdtree_new(numverts);
+	for(i = 0; i < numfaces; i++)
+	{
+		VECCOPY(co1, mvert[mface[i].v1].co);
+		Mat4MulVecfl(mat, co1);
+		set_min_max(min, max, co1);
+		VECCOPY(co2, mvert[mface[i].v2].co);
+		Mat4MulVecfl(mat, co2);
+		set_min_max(min, max, co2);
+		VECCOPY(co3, mvert[mface[i].v3].co);
+		Mat4MulVecfl(mat, co3);
+		set_min_max(min, max, co3);
+		
+		// calc triangle center
+		VECCOPY(co, co1);
+		VECADD(co, co, co2);
+		VECADD(co, co, co3);
+		if(mface[i].v4)
+		{
+			VECCOPY(co4, mvert[mface[i].v4].co);
+			Mat4MulVecfl(mat, co4);
+			set_min_max(min, max, co4);
+			VECADD(co, co, co4);
+			VecMulf(co, 0.25);
+		}
+		else
+			VecMulf(co, 1.0 / 3.0);
+		
+		if(mface[i].v4)
+			CalcNormFloat4(mvert[mface[i].v1].co, mvert[mface[i].v2].co, mvert[mface[i].v3].co, mvert[mface[i].v4].co, &normals[i*3]);
+		else
+			CalcNormFloat(mvert[mface[i].v1].co, mvert[mface[i].v2].co, mvert[mface[i].v3].co, &normals[i*3]);
+		
+		BLI_kdtree_insert(tree, i, co, NULL);
+	}
+	BLI_kdtree_balance(tree);
+	/////////////////////////////////////////////////
+	
+	// calculate slice height + width
+	maxdir = max[0] - min[0];
+	maxdir = MAX2(max[1]-min[1], maxdir);
+	maxdir = MAX2(max[2]-min[2], maxdir);
+	slice = maxdir / (float)maxres;
+	resx = MIN2(maxres, ceil((max[0] - min[0]) / slice));
+	resy = MIN2(maxres, ceil((max[1] - min[1]) / slice));
+	resz = MIN2(maxres, ceil((max[2] - min[2]) / slice));
+	
+	// adjust max
+	max[0] = min[0] + slice * resx;
+	max[1] = min[1] + slice * resy;
+	max[2] = min[2] + slice * resz;
+	
+	if(sphmd->sim_parms->co)
+		MEM_freeN(sphmd->sim_parms->co);
+	
+	cos = sphmd->sim_parms->co = MEM_callocN(sizeof(float)*3*resx*resy*resz, "sph_cos");
+	// r = calloc(1, sizeof(float)*resx*resy*resz);
+	maxpart = sizeof(float)*3*resx*resy*resz;
+	
+// #pragma omp parallel for private(i,j,k) schedule(static)
+	for(i = 0; i < resx; i ++)
+	{
+		for(j = 0; j < resy; j++)
+		{
+			for(k = 0; k < resz; k++)
+			{
+				KDTreeNearest nearest;
+				float tco[3];
+				tco[0] = min[0] + slice * i + slice * 0.5;
+				tco[1] = min[1] + slice * j + slice * 0.5;
+				tco[2] = min[2] + slice * k + slice * 0.5;
+				
+				index = BLI_kdtree_find_nearest(tree, tco, NULL, &nearest);
+				
+				if(index != -1)
+				{
+					float t[3];
+					float sgn;
+					
+					VECSUB(t, tco, nearest.co);
+					sgn = INPR(t, &normals[nearest.index*3]);
+					
+					if(sgn < 0.0)
+						sgn = -1.0;
+					else 
+						sgn = 1.0;
+					
+					dist[(i*resy*resz)+(j*resz)+k] = sgn * nearest.dist;
+					
+					if((int)sgn < 0)
+					{	
+						// create particle if inside object
+						VECCOPY(&cos[totpart*3], tco);
+						
+						totpart++;	
+					}
+					
+				}
+				else
+				{
+					printf("Error: no nearest point!\n");
+				}
+			}
+		}
+	}
+	printf("maxpart: %d, totpart: %d\n", maxpart, totpart);
+	
+	sphmd->sim_parms->numpart = totpart;
+	MEM_freeN(dist);
+	MEM_freeN(normals);
+	BLI_kdtree_free(tree);
+	return totpart;
+}
+
+/* add constraints, inflows, fluid, ... */
+int sph_init_all (SphModifierData *sphmd, DerivedMesh *dm, Object *ob)
+{
+	Base *base=NULL;
+	Object *fobject = NULL;
+	SphModifierData *fsphmd = NULL;
+	int fluids = 0; // only one fluid object possible
+	DerivedMesh *fdm = NULL;
+	
+	sphmd->sim_parms->flags |= SPH_SIMSETTINGS_FLAG_INIT;
+	CDDM_calc_normals ( dm );
+	
+	/* create C++ object */
+	// sph_init_cpp(sphmd);
+	sphmd->sim_parms->numpart = calc_distance_field(ob, dm, sphmd, ob->obmat);
+	return 1;
+	
+	/* create fluid domain */
+	// sph_set_domain(sphmd, dm, ob->obmat);
+	
+	// check for constraints, fluid, etc. but ignore domains
+	for ( base = G.scene->base.first; base; base = base->next )
+	{
+		fobject = base->object;
+		fsphmd = ( SphModifierData * ) modifiers_findByType ( fobject, eModifierType_Sph );
+		
+		// TODO I could check for linked groups, too 
+		if ( !fsphmd )
+		{
+			// TODO
+		}
+		else
+		{	
+			if(fsphmd == sphmd)
+				continue;
+			
+			if(fsphmd->sim_parms)
+			{	
+				// check for fluid
+				if(((short)fsphmd->sim_parms->flags == SPH_SIMSETTINGS_FLAG_FLUID)  && (!fluids))
+				{
+					// create fluids
+					// particles have to be created AFTER constraints
+					// TODO: no particles = crash
+					
+					// get derivedmesh from object
+					fdm = mesh_get_derived_final(fobject, CD_MASK_BAREMESH);
+					
+					// create fluid object
+					/*
+					if(!sph_add_particles(sphmd, fdm, fobject->obmat, fsphmd->sim_parms->resolution))
+					{
+						fluids--;
+					}
+					*/
+					if(fdm)
+						fdm->release(fdm);
+					
+					fluids++;
+				}
+				else if((short)fsphmd->sim_parms->flags & SPH_SIMSETTINGS_FLAG_OBSTACLE)
+				{
+					
+				}
+			}
+		}
+	}
+	
+	if(!fluids)
+		return 0;
+		
+	return 1;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h
new file mode 100644
index 00000000000..51f87b26aaf
--- /dev/null
+++ b/source/blender/blenlib/BLI_kdopbvh.h
@@ -0,0 +1,63 @@
+/**
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * 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 LICENSE BLOCK *****
+ */
+
+
+#ifndef BLI_KDOPBVH_H
+#define BLI_KDOPBVH_H
+
+#include <float.h>
+
+struct BVHTree;
+typedef struct BVHTree BVHTree;
+
+typedef struct BVHTreeOverlap {
+	int indexA;
+	int indexB;
+} BVHTreeOverlap;
+
+BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis);
+void BLI_bvhtree_free(BVHTree *tree);
+
+/* construct: first insert points, then call balance */
+int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints);
+void BLI_bvhtree_balance(BVHTree *tree);
+
+/* update: first update points/nodes, then call update_tree to refit the bounding volumes */
+int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints);
+void BLI_bvhtree_update_tree(BVHTree *tree);
+
+/* collision/overlap: check two trees if they overlap, alloc's *overlap with length of the int return value */
+BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result);
+
+float BLI_bvhtree_getepsilon(BVHTree *tree);
+
+#endif // BLI_KDOPBVH_H
+
+
+
+
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
new file mode 100644
index 00000000000..8be52854a7b
--- /dev/null
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -0,0 +1,786 @@
+/*  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): Daniel Genrich
+*
+* ***** END GPL/BL DUAL LICENSE BLOCK *****
+*/
+
+#include "math.h"
+#include <stdio.h>
+#include <stdlib.h> 
+#include <string.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_utildefines.h"
+
+#include "BLI_kdopbvh.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+typedef struct BVHNode
+{
+	struct BVHNode *children[8]; // max 8 children
+	struct BVHNode *parent; // needed for bottom - top update
+	float bv[26]; // Bounding volume of all nodes, max 13 axis
+	int index; /* face, edge, vertex index */
+	char totnode; // how many nodes are used, used for speedup
+	char traversed;  // how many nodes already traversed until this level?
+	char main_axis;
+} BVHNode;
+
+struct BVHTree
+{
+	BVHNode **nodes;
+	BVHNode *nodearray; /* pre-alloc branch nodes */
+	int	*orig_index; /* mapping for orig_index to node_index */
+	float 	epsilon; /* epslion is used for inflation of the k-dop	   */
+	int 	totleaf; // leafs
+	int 	totbranch;
+	char 	tree_type; // type of tree (4 => quadtree)
+	char 	axis; // kdop type (6 => OBB, 7 => AABB, ...)
+	char 	start_axis, stop_axis; // KDOP_AXES array indices according to axis
+};
+
+typedef struct BVHOverlapData 
+{  
+	BVHTree *tree1, *tree2; 
+	BVHTreeOverlap *overlap; 
+	int i, max_overlap; /* i is number of overlaps */
+} BVHOverlapData;
+////////////////////////////////////////
+
+
+////////////////////////////////////////////////////////////////////////
+// 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}, {0, 1.0, 0}, {0, 0, 1.0}, {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0},
+{1.0, -1.0, -1.0}, {1.0, 1.0, 0}, {1.0, 0, 1.0}, {0, 1.0, 1.0}, {1.0, -1.0, 0}, {1.0, 0, -1.0},
+{0, 1.0, -1.0}
+};
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+// 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
+*/
+static void bvh_exchange(BVHNode **a, int i, int j)
+{
+	BVHNode *t=a[i];
+	a[i]=a[j];
+	a[j]=t;
+}
+static int floor_lg(int a)
+{
+	return (int)(floor(log(a)/log(2)));
+}
+
+/*
+* Insertion sort algorithm
+*/
+static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
+{
+	int i,j;
+	BVHNode *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(BVHNode **a, int lo, int hi, BVHNode * 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(BVHNode **a, int i, int n, int lo, int axis)
+{
+	BVHNode * 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(BVHNode **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 BVHNode *bvh_medianof3(BVHNode **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 (BVHNode **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;
+	}
+}
+
+static void sort(BVHNode **a0, int begin, int end, int axis)
+{
+	if (begin < end)
+	{
+		BVHNode **a=a0;
+		bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
+		bvh_insertionsort(a, begin, end, axis);
+	}
+}
+void sort_along_axis(BVHTree *tree, int start, int end, int axis)
+{
+	sort(tree->nodes, start, end, axis);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void BLI_bvhtree_free(BVHTree *tree)
+{	
+	if(tree)
+	{
+		MEM_freeN(tree->nodes);
+		MEM_freeN(tree->nodearray);
+		MEM_freeN(tree->orig_index);
+		MEM_freeN(tree);
+	}
+}
+
+BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
+{
+	BVHTree *tree;
+	int numbranches=0, i;
+	
+	// only support up to octree
+	if(tree_type > 8)
+		return NULL;
+
+	tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
+	
+	if(tree)
+	{
+		// calculate max number of branches, our bvh kdop is "almost perfect"
+		for(i = 1; i <= (int)ceil((float)((float)log(maxsize)/(float)log(tree_type))); i++)
+			numbranches += (pow(tree_type, i) / tree_type);
+		
+		tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*(numbranches+maxsize + tree_type), "BVHNodes");
+		
+		if(!tree->nodes)
+		{
+			MEM_freeN(tree);
+			return NULL;
+		}
+		
+		tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)*(numbranches+maxsize + tree_type), "BVHNodeArray");
+		
+		if(!tree->nodearray)
+		{
+			MEM_freeN(tree);
+			MEM_freeN(tree->nodes);
+			return NULL;
+		}
+		
+		tree->orig_index = (int *)MEM_callocN(sizeof(int)*(numbranches+maxsize + tree_type), "BVHIndexArray");
+		
+		if(!tree->orig_index)
+		{
+			MEM_freeN(tree);
+			MEM_freeN(tree->nodes);
+			MEM_freeN(tree->nodearray);
+			return NULL;
+		}
+		
+		tree->epsilon = epsilon;
+		tree->tree_type = tree_type; 
+		tree->axis = axis;
+		
+		if(axis == 26)
+		{
+			tree->start_axis = 0;
+			tree->stop_axis = 13;
+		}
+		else if(axis == 18)
+		{
+			tree->start_axis = 7;
+			tree->stop_axis = 13;
+		}
+		else if(axis == 14)
+		{
+			tree->start_axis = 0;
+			tree->stop_axis = 7;
+		}
+		else if(axis == 8) // AABB
+		{
+			tree->start_axis = 0;
+			tree->stop_axis = 4;
+		}
+		else if(axis == 6) // OBB
+		{
+			tree->start_axis = 0;
+			tree->stop_axis = 3;
+		}
+		else
+		{
+			BLI_bvhtree_free(tree);
+			return NULL;
+		}
+	}
+
+	return tree;
+}
+
+static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoints, int moving)
+{
+	float newminmax;
+	int i, k;
+	
+	// don't init boudings for the moving case
+	if(!moving)
+	{
+		for (i = tree->start_axis; i < tree->stop_axis; i++)
+		{
+			node->bv[2*i] = FLT_MAX;
+			node->bv[2*i + 1] = -FLT_MAX;
+		}
+	}
+	
+	for(k = 0; k < numpoints; k++)
+	{
+		// for all Axes.
+		for (i = tree->start_axis; i < tree->stop_axis; i++)
+		{
+			newminmax = INPR(&co[k * 3], KDOP_AXES[i]);
+			if (newminmax < node->bv[2 * i])
+				node->bv[2 * i] = newminmax;
+			if (newminmax > node->bv[(2 * i) + 1])
+				node->bv[(2 * i) + 1] = newminmax;
+		}
+	}
+}
+
+// depends on the fact that the BVH's for each face is already build
+static void refit_kdop_hull(BVHTree *tree, int start, int end, float *bv)
+{
+	float newmin,newmax;
+	int i, j;
+	
+	for (i = tree->start_axis; i < tree->stop_axis; i++)
+	{
+		bv[2*i] = FLT_MAX;
+		bv[2*i + 1] = -FLT_MAX;
+	}
+
+	for (j = start; j < end; j++)
+	{
+                // for all Axes.
+		for (i = tree->start_axis; i < tree->stop_axis; i++)
+		{
+			newmin = tree->nodes[j]->bv[(2 * i)];   
+			if ((newmin < bv[(2 * i)]))
+				bv[(2 * i)] = newmin;
+ 
+			newmax = tree->nodes[j]->bv[(2 * i) + 1];
+			if ((newmax > bv[(2 * i) + 1]))
+				bv[(2 * i) + 1] = newmax;
+		}
+	}
+}
+
+int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints)
+{
+	BVHNode *node= NULL;
+	int i;
+	
+	// insert should only possible as long as tree->totbranch is 0
+	if(tree->totbranch > 0)
+		return 0;
+	
+	if(tree->totleaf+1 >= MEM_allocN_len(tree->nodes))
+		return 0;
+	
+	// TODO check if have enough nodes in array
+	
+	node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
+	tree->totleaf++;
+	
+	create_kdop_hull(tree, node, co, numpoints, 0);
+	
+	// inflate the bv with some epsilon
+	for (i = tree->start_axis; i < tree->stop_axis; i++)
+	{
+		node->bv[(2 * i)] -= tree->epsilon; // minimum 
+		node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
+	}
+
+	node->index= index;
+	
+	return 1;
+}
+
+// only supports x,y,z axis in the moment
+// but we should use a plain and simple function here for speed sake
+static char get_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
+	}
+}
+
+static void bvh_div_nodes(BVHTree *tree, BVHNode *node, int start, int end, char lastaxis)
+{
+	char laxis;
+	int i, tend;
+	BVHNode *tnode;
+	int slice = (end-start+tree->tree_type-1)/tree->tree_type;	//division rounded up
+	
+	// Determine which axis to split along
+	laxis = get_largest_axis(node->bv);
+
+	// Sort along longest axis
+	if(laxis!=lastaxis)
+		sort_along_axis(tree, start, end, laxis);
+	
+	// split nodes along longest axis
+	for (i=0; start < end; start += slice, i++) //i counts the current child
+	{	
+		tend = start + slice;
+		
+		if(tend > end) tend = end;
+		
+		if(tend-start == 1)	// ok, we have 1 left for this node
+		{
+			node->children[i] = tree->nodes[start];
+			node->children[i]->parent = node;
+		}
+		else
+		{
+			tnode = node->children[i] = tree->nodes[tree->totleaf  + tree->totbranch] = &(tree->nodearray[tree->totbranch + tree->totleaf]);
+			tree->totbranch++;
+			tnode->parent = node;
+
+			refit_kdop_hull(tree, start, tend, tnode->bv);
+			bvh_div_nodes(tree, tnode, start, tend, laxis);
+		}
+		node->totnode++;
+	}
+	
+	return;
+}
+
+static void verify_tree(BVHTree *tree)
+{
+	int i, j, check = 0;
+	
+	// check the pointer list
+	for(i = 0; i < tree->totleaf; i++)
+	{
+		if(tree->nodes[i]->parent == NULL)
+			printf("Leaf has no parent: %d\n", i);
+		else
+		{
+			for(j = 0; j < tree->tree_type; j++)
+			{
+				if(tree->nodes[i]->parent->children[j] == tree->nodes[i])
+					check = 1;
+			}
+			if(!check)
+			{
+				printf("Parent child relationship doesn't match: %d\n", i);
+			}
+			check = 0;
+		}
+	}
+	
+	// check the leaf list
+	for(i = 0; i < tree->totleaf; i++)
+	{
+		if(tree->nodearray[i].parent == NULL)
+			printf("Leaf has no parent: %d\n", i);
+		else
+		{
+			for(j = 0; j < tree->tree_type; j++)
+			{
+				if(tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
+					check = 1;
+			}
+			if(!check)
+			{
+				printf("Parent child relationship doesn't match: %d\n", i);
+			}
+			check = 0;
+		}
+	}
+	
+	printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
+}
+	
+void BLI_bvhtree_balance(BVHTree *tree)
+{
+	BVHNode *node;
+	int i;
+	
+	if(tree->totleaf == 0)
+		return;
+	
+	// create root node
+	node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
+	tree->totbranch++;
+	
+	// refit root bvh node
+	refit_kdop_hull(tree, 0, tree->totleaf, tree->nodes[tree->totleaf]->bv);
+	// create + balance tree
+	bvh_div_nodes(tree, tree->nodes[tree->totleaf], 0, tree->totleaf, 0);
+	
+	// put indices into array for O(1) access
+	for(i = 0; i < tree->totleaf; i++)
+	{
+		tree->orig_index[tree->nodes[i]->index] = i;
+	}
+	
+	verify_tree(tree);
+}
+
+// overlap - is it possbile for 2 bv's to collide ?
+static int tree_overlap(float *bv1, float *bv2, int start_axis, int stop_axis)
+{
+	float *bv1_end = bv1 + (stop_axis<<1);
+		
+	bv1 += start_axis<<1;
+	bv2 += start_axis<<1;
+	
+	// test all axis if min + max overlap
+	for (; bv1 != bv1_end; bv1+=2, bv2+=2)
+	{
+		if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
+			return 0;
+	}
+	
+	return 1;
+}
+
+static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
+{
+	int j;
+	
+	if(tree_overlap(node1->bv, node2->bv, MIN2(data->tree1->start_axis, data->tree2->start_axis), MIN2(data->tree1->stop_axis, data->tree2->stop_axis)))
+	{
+		// check if node1 is a leaf
+		if(node1->index)
+		{
+			// check if node2 is a leaf
+			if(node2->index)
+			{
+				
+				if(node1 == node2)
+				{
+					return;
+				}
+					
+				if(data->i >= data->max_overlap)
+				{	
+					// try to make alloc'ed memory bigger
+					data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
+					
+					if(!data->overlap)
+					{
+						printf("Out of Memory in traverse\n");
+						return;
+					}
+					data->max_overlap *= 2;
+				}
+				
+				// both leafs, insert overlap!
+				data->overlap[data->i].indexA = node1->index;
+				data->overlap[data->i].indexB = node2->index;
+
+				data->i++;
+			}
+			else
+			{
+				for(j = 0; j < data->tree2->tree_type; j++)
+				{
+					if(node2->children[j])
+						traverse(data, node1, node2->children[j]);
+				}
+			}
+		}
+		else
+		{
+			
+			for(j = 0; j < data->tree2->tree_type; j++)
+			{
+				if(node1->children[j])
+					traverse(data, node1->children[j], node2);
+			}
+		}
+	}
+	return;
+}
+
+BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result)
+{
+	int j, total = 0;
+	BVHTreeOverlap *overlap = NULL, *to = NULL;
+	BVHOverlapData *data[tree1->tree_type];
+	
+	// check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
+	if((tree1->axis != tree2->axis) && ((tree1->axis == 14) || tree2->axis == 14))
+		return 0;
+	
+	// fast check root nodes for collision before doing big splitting + traversal
+	if(!tree_overlap(tree1->nodes[tree1->totleaf]->bv, tree2->nodes[tree2->totleaf]->bv, MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis)))
+		return 0;
+	
+	for(j = 0; j < tree1->tree_type; j++)
+	{
+		data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
+		
+		// init BVHOverlapData
+		data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
+		data[j]->tree1 = tree1;
+		data[j]->tree2 = tree2;
+		data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
+		data[j]->i = 0;
+	}
+
+#pragma omp parallel for private(j) schedule(static)
+	for(j = 0; j < tree1->tree_type; j++)
+	{
+		traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
+	}
+	
+	for(j = 0; j < tree1->tree_type; j++)
+		total += data[j]->i;
+	
+	to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
+	
+	for(j = 0; j < tree1->tree_type; j++)
+	{
+		memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
+		to+=data[j]->i;
+	}
+	
+	for(j = 0; j < tree1->tree_type; j++)
+	{
+		free(data[j]->overlap);
+		MEM_freeN(data[j]);
+	}
+	
+	(*result) = total;
+	return overlap;
+}
+
+
+// bottom up update of bvh tree:
+// join the 4 children here
+static void node_join(BVHTree *tree, BVHNode *node)
+{
+	int i, j;
+	
+	for (i = tree->start_axis; i < tree->stop_axis; i++)
+	{
+		node->bv[2*i] = FLT_MAX;
+		node->bv[2*i + 1] = -FLT_MAX;
+	}
+	
+	for (i = 0; i < tree->tree_type; i++)
+	{
+		if (node->children[i]) 
+		{
+			for (j = tree->start_axis; j < tree->stop_axis; j++)
+			{
+				// update minimum 
+				if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
+					node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
+				
+				// update maximum 
+				if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
+					node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
+			}
+		}
+		else
+			break;
+	}
+}
+
+// call before BLI_bvhtree_update_tree()
+int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints)
+{
+	BVHNode *node= NULL;
+	int i = 0;
+	
+	// check if index exists
+	if(index > tree->totleaf)
+		return 0;
+	
+	node = tree->nodes[tree->orig_index[index]];
+	
+	create_kdop_hull(tree, node, co, numpoints, 0);
+	
+	if(co_moving)
+		create_kdop_hull(tree, node, co_moving, numpoints, 1);
+	
+	// inflate the bv with some epsilon
+	for (i = tree->start_axis; i < tree->stop_axis; i++)
+	{
+		node->bv[(2 * i)] -= tree->epsilon; // minimum 
+		node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
+	}
+	
+	return 1;
+}
+
+// call BLI_bvhtree_update_node() first for every node/point/triangle
+void BLI_bvhtree_update_tree(BVHTree *tree)
+{
+	BVHNode *leaf, *parent;
+	
+	// reset tree traversing flag
+	for (leaf = tree->nodearray + tree->totleaf; leaf != tree->nodearray + tree->totleaf + tree->totbranch; leaf++)
+		leaf->traversed = 0;
+	
+	for (leaf = tree->nodearray; leaf != tree->nodearray + tree->totleaf; leaf++)
+	{
+		for (parent = leaf->parent; parent; parent = parent->parent)
+		{
+			parent->traversed++;	// we tried to go up in hierarchy 
+			if (parent->traversed < parent->totnode) 
+				break;	// we do not need to check further 
+			else 
+				node_join(tree, parent);
+		}
+	}
+}
+
+float BLI_bvhtree_getepsilon(BVHTree *tree)
+{
+	return tree->epsilon;
+}
diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c
index 24606bd822f..e3649a0cc43 100644
--- a/source/blender/blenloader/intern/readfile.c
+++ b/source/blender/blenloader/intern/readfile.c
@@ -3073,6 +3073,7 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 			
 			clmd->sim_parms= newdataadr(fd, clmd->sim_parms);
 			clmd->coll_parms= newdataadr(fd, clmd->coll_parms);
+
 			clmd->point_cache= newdataadr(fd, clmd->point_cache);
 
 			if(clmd->point_cache)
@@ -3082,22 +3083,16 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 				if(clmd->sim_parms->presets > 10)
 					clmd->sim_parms->presets = 0;
 			}
+			else
+			{
+				/* Collision modifier without parameters?? */
+				
+			}
 			
 		}
 		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;
@@ -3106,9 +3101,8 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
 			collmd->current_v = NULL;
 			collmd->time = -1;
 			collmd->numverts = 0;
-			collmd->bvh = NULL;
+			collmd->bvhtree = NULL;
 			collmd->mfaces = NULL;
-			
 		}
 		else if (md->type==eModifierType_Hook) {
 			HookModifierData *hmd = (HookModifierData*) md;
diff --git a/source/blender/include/butspace.h b/source/blender/include/butspace.h
index 242965a820c..c6482290c86 100644
--- a/source/blender/include/butspace.h
+++ b/source/blender/include/butspace.h
@@ -297,6 +297,10 @@ void curvemap_buttons(struct uiBlock *block, struct CurveMapping *cumap, char la
 /* Cloth sim button defines */
 #define B_CLOTH_CHANGEPREROLL	1480
 
+/* SPH sim button defines */
+#define B_SPH_BAKE 		1490
+
+
 /* *********************** */
 #define B_WORLDBUTS		1600
 
diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h
index b7b43817474..b2d6e8e5756 100644
--- a/source/blender/makesdna/DNA_modifier_types.h
+++ b/source/blender/makesdna/DNA_modifier_types.h
@@ -35,6 +35,7 @@ typedef enum ModifierType {
 	eModifierType_Cloth,
 	eModifierType_Collision,
 	eModifierType_Bevel,
+	eModifierType_Sph,
 	NUM_MODIFIER_TYPES
 } ModifierType;
 
@@ -374,6 +375,14 @@ typedef struct ClothModifierData {
    struct PointCache *point_cache;	/* definition is in DNA_object_force.h */
 } ClothModifierData;
 
+typedef struct SphModifierData {
+	ModifierData		modifier;
+	
+	struct SPH *sph; /* pointer to cpp sph sim */
+	struct SphSimSettings *sim_parms; /* definition is in DNA_sph_types.h */
+	struct SphCollSettings *coll_parms; /* definition is in DNA_sph_types.h */
+} SphModifierData;
+
 typedef struct CollisionModifierData {
 	ModifierData	modifier;
 	
@@ -390,7 +399,7 @@ typedef struct CollisionModifierData {
 	unsigned int numfaces;
 	int pad;
 	float time;		/* cfra time of modifier */
-	struct BVH *bvh;	/* bounding volume hierarchy for this cloth object */
+	struct BVHTree *bvhtree; /* bounding volume hierarchy for this cloth object */
 } CollisionModifierData;
 
 typedef enum {
diff --git a/source/blender/makesdna/DNA_sph_types.h b/source/blender/makesdna/DNA_sph_types.h
new file mode 100644
index 00000000000..add3053f612
--- /dev/null
+++ b/source/blender/makesdna/DNA_sph_types.h
@@ -0,0 +1,102 @@
+/**
+ * $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_SPH_TYPES_H
+#define DNA_SPH_TYPES_H
+
+// string scene, constraint, backup: missing
+//	string _file;
+//	string pov;
+//	string rendered;
+typedef struct SphSimSettings
+{
+	int flags; // see pw_extern.h
+	float timestep; 
+	float viscosity;
+	float incompressibility; /* how incompressible is the fluid? */
+	float surfacetension;
+	float density;
+	float gravity[3]; /* gravity on the domain */
+	float samplingdistance;
+	float smoothinglength;
+
+	float controlviscosity;
+	int dumpimageevery;
+	int computesurfaceevery;
+	int fastmarchingevery;
+	int dumppovrayevery;
+	
+	float totaltime;
+	
+	float tangentialfriction; /* constraint tangential friction */
+	float normalfriction; /* constraint normal friction */
+
+	float rotation_angle;
+	float rotation_axis[3];
+	float rotation_center[3];
+	float scenelowerbound[3];
+	float sceneupperbound[3];
+	
+	int initiallevel;
+	
+	float alpha;
+	float beta;
+	float gamma;
+	
+	/* needed for better direct resolution input for constraints, fluids,... */
+	int resolution; /* can also be calculated by (Max-Min) / samplingdistance */
+	int pad;
+	float *verts;
+	float *normals;
+	int *tris;
+	unsigned int numverts;
+	unsigned int numtris;
+	float *co; /* particle positions */
+	float *r; /* particle radius */
+	long numpart;
+	int pad2;
+}
+SphSimSettings;
+
+typedef struct SphCollSettings
+{
+	float	epsilon;		/* min distance for collisions.		*/
+	float	self_friction;		/* Fiction/damping with self contact.		 	*/
+	float	friction;		/* Friction/damping applied on contact with other object.*/
+	short	self_loop_count;	/* How many iterations for the selfcollision loop	*/
+	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 	selfepsilon; 		/* for selfcollision */
+}
+SphCollSettings;
+
+#endif // DNA_SPH_TYPES_H
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index b02e89233e5..815a9c3869e 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -985,7 +985,7 @@ static uiBlock *modifiers_add_menu(void *ob_v)
 		/* Only allow adding through appropriate other interfaces */
 		if(ELEM3(i, eModifierType_Softbody, eModifierType_Hook, eModifierType_ParticleSystem)) continue;
 		
-		if(ELEM(i, eModifierType_Cloth, eModifierType_Collision)) continue;
+		if(ELEM3(i, eModifierType_Cloth, eModifierType_Collision, eModifierType_Sph)) continue;
 
 		if((mti->flags&eModifierTypeFlag_AcceptsCVs) ||
 		   (ob->type==OB_MESH && (mti->flags&eModifierTypeFlag_AcceptsMesh))) {
@@ -1807,6 +1807,8 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			height = 31;
 		} else if (md->type==eModifierType_Collision) {
 			height = 31;
+		} else if (md->type==eModifierType_Sph) {
+			height = 31;
 		} else if (md->type==eModifierType_Boolean) {
 			height = 48;
 		} else if (md->type==eModifierType_Array) {
@@ -1837,7 +1839,7 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 				uiButSetFunc(but, modifiers_applyModifier, ob, md);
 			}
 			
-			if (md->type!=eModifierType_Softbody && md->type!=eModifierType_ParticleSystem && (md->type!=eModifierType_Cloth)) {
+			if (md->type!=eModifierType_Softbody && md->type!=eModifierType_ParticleSystem && (md->type!=eModifierType_Cloth) && (md->type!=eModifierType_Sph)) {
 				but = uiDefBut(block, BUT, B_MODIFIER_RECALC, "Copy",	lx,(cy-=19),60,19, 0, 0, 0, 0, 0, "Duplicate the current modifier at the same position in the stack");
 				uiButSetFunc(but, modifiers_copyModifier, ob, md);
 			}
@@ -2229,6 +2231,8 @@ static void draw_modifier(uiBlock *block, Object *ob, ModifierData *md, int *xco
 			uiDefBut(block, LABEL, 1, "See Soft Body panel.",	lx, (cy-=19), buttonWidth,19, NULL, 0.0, 0.0, 0, 0, "");
 		} else if (md->type==eModifierType_Cloth) {
 			uiDefBut(block, LABEL, 1, "See Cloth panel.",	lx, (cy-=19), buttonWidth,19, NULL, 0.0, 0.0, 0, 0, "");
+		} else if (md->type==eModifierType_Sph) {
+			uiDefBut(block, LABEL, 1, "See Sph panel.",	lx, (cy-=19), buttonWidth,19, NULL, 0.0, 0.0, 0, 0, "");
 		} else if (md->type==eModifierType_Collision) {
 			uiDefBut(block, LABEL, 1, "See Collision panel.",	lx, (cy-=19), buttonWidth,19, NULL, 0.0, 0.0, 0, 0, "");
 		} else if (md->type==eModifierType_Boolean) {
diff --git a/source/blender/src/buttons_object.c b/source/blender/src/buttons_object.c
index 9a749b6fbd9..4dcc51e03e3 100644
--- a/source/blender/src/buttons_object.c
+++ b/source/blender/src/buttons_object.c
@@ -52,6 +52,8 @@
 #include "BKE_utildefines.h"
 #include "BKE_particle.h"
 #include "BKE_pointcache.h"
+#include "BKE_sph.h"
+
 
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
@@ -107,6 +109,7 @@
 #include "DNA_particle_types.h"
 #include "DNA_radio_types.h"
 #include "DNA_screen_types.h"
+#include "DNA_sph_types.h"
 #include "DNA_sound_types.h"
 #include "DNA_texture_types.h"
 #include "DNA_userdef_types.h"
@@ -135,6 +138,7 @@
 #include "BKE_object.h"
 #include "BKE_particle.h"
 #include "BKE_sound.h"
+#include "BKE_sph.h"
 #include "BKE_texture.h"
 #include "BKE_utildefines.h"
 #include "BKE_DerivedMesh.h"
@@ -2476,6 +2480,27 @@ void do_object_panels(unsigned short event)
 		allqueue(REDRAWVIEW3D, 0);
 	}
 	break;
+	case B_SPH_BAKE:
+	{
+		SphModifierData *sphmd = (SphModifierData *)modifiers_findByType(ob, eModifierType_Sph);
+		int i = 0;
+		
+		if(sphmd && sphmd->sim_parms)
+		{
+			sphmd->sim_parms->flags |= SPH_SIMSETTINGS_FLAG_BAKING;
+			
+			// call baking function
+			for(i = 0; i < 1; i++)
+			{
+				CFRA++;
+				update_for_newframe();
+			}
+			
+			sphmd->sim_parms->flags &= ~SPH_SIMSETTINGS_FLAG_BAKING;
+			
+		}
+	}
+	break;
 	}
 
 }
@@ -5650,6 +5675,116 @@ static void object_panel_cloth_III(Object *ob)
 	
 }
 
+static void object_sph__enabletoggle(void *ob_v, void *arg2)
+{
+	Object *ob = ob_v;
+	ModifierData *md = modifiers_findByType(ob, eModifierType_Sph);
+
+	if (!md) {
+		// create particle modifier
+		ParticleSettings *part = psys_new_settings("PSys", G.main);
+		ParticleSystem *psys = MEM_callocN(sizeof(ParticleSystem), "particle_system");
+		ParticleSystemModifierData *psmd;
+		
+		md = modifier_new(eModifierType_Sph);
+		BLI_addtail(&ob->modifiers, md);
+
+		part->type = PART_FLUID;
+		psys->part = part;
+		psys->pointcache = BKE_ptcache_add();
+		psys->flag |= PSYS_ENABLED;
+		BLI_addtail(&ob->particlesystem,psys);
+		md= modifier_new(eModifierType_ParticleSystem);
+		sprintf(md->name, "SphParticleSystem" );
+		psmd= (ParticleSystemModifierData*) md;
+		psmd->psys=psys;
+		BLI_addtail(&ob->modifiers, md);
+		
+		DAG_object_flush_update(G.scene, ob, OB_RECALC_DATA);
+		allqueue(REDRAWBUTSEDIT, 0);
+		allqueue(REDRAWVIEW3D, 0);
+	}
+	else {
+		Object *ob = ob_v;
+		ModifierData *md = modifiers_findByType(ob, eModifierType_Sph);
+	
+		if (!md)
+			return;
+
+		BLI_remlink(&ob->modifiers, md);
+
+		modifier_free(md);
+
+		BIF_undo_push("Del modifier");
+		
+		// ob->softflag |= OB_SB_RESET;
+		allqueue(REDRAWBUTSEDIT, 0);
+		allqueue(REDRAWVIEW3D, 0);
+		allqueue(REDRAWIMAGE, 0);
+		allqueue(REDRAWOOPS, 0);
+		DAG_object_flush_update(G.scene, ob, OB_RECALC_DATA);
+		object_handle_update(ob);
+		countall();
+	}
+}
+
+static void object_panel_sph(Object *ob)
+{
+	uiBlock *block=NULL;
+	uiBut *but=NULL;
+	static int val;
+	SphModifierData *sphmd = (SphModifierData *)modifiers_findByType(ob, eModifierType_Sph);
+	
+	block= uiNewBlock(&curarea->uiblocks, "object_sph", UI_EMBOSS, UI_HELV, curarea->win);
+	if(uiNewPanel(curarea, block, "Sph ", "Physics", 640, 0, 318, 204)==0) return;
+	uiSetButLock(object_data_is_libdata(ob), ERROR_LIBDATA_MESSAGE);
+	
+	val = (sphmd ? 1:0);
+	
+
+	but = uiDefButI(block, TOG, REDRAWBUTSOBJECT, "Sph",	10,200,130,20, &val, 0, 0, 0, 0, "Sets object to become Sph");
+	uiButSetFunc(but, object_sph__enabletoggle, ob, NULL);
+
+
+	uiDefBut(block, LABEL, 0, "",10,10,300,0, NULL, 0.0, 0, 0, 0, ""); /* tell UI we go to 10,10*/
+	
+	if(sphmd)
+	{
+		if(sphmd->sim_parms && (sphmd->sim_parms->flags & SPH_SIMSETTINGS_FLAG_DOMAIN))
+		{
+			uiDefBut(block, BUT, B_SPH_BAKE, "BAKE",10, 180,300,20, NULL, 0.0, 0.0, 10, 0, "Perform simulation and output and surface&preview meshes for each frame.");
+		}
+		
+		uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Domain", 10, 160, 100, 20, (short *)&sphmd->sim_parms->flags, 15.0, SPH_SIMSETTINGS_FLAG_DOMAIN,   0.0, 0.0, " ");
+		uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Fluid", 110, 160, 100, 20, (short *)&sphmd->sim_parms->flags, 15.0, SPH_SIMSETTINGS_FLAG_FLUID,   0.0, 0.0, " ");
+		uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Obstacle", 210, 160, 100, 20, (short *)&sphmd->sim_parms->flags, 15.0, SPH_SIMSETTINGS_FLAG_OBSTACLE,   0.0, 0.0, " ");
+		
+		if(sphmd->sim_parms && (sphmd->sim_parms->flags & SPH_SIMSETTINGS_FLAG_DOMAIN))
+		{
+			uiDefButBitI(block, TOG, SPH_SIMSETTINGS_FLAG_GHOSTS, REDRAWBUTSOBJECT, "Ghosts",10,140,100,20, &sphmd->sim_parms->flags, 0, 0, 0, 0, " ");
+			uiDefButBitI(block, TOG, SPH_SIMSETTINGS_FLAG_MULTIRES, REDRAWBUTSOBJECT, "Multires",110,140,100,20, &sphmd->sim_parms->flags, 0, 0, 0, 0, " ");
+			uiDefButBitI(block, TOG, SPH_SIMSETTINGS_FLAG_VORTICITY, REDRAWBUTSOBJECT, "Vorticity",210,140,100,20, &sphmd->sim_parms->flags, 0, 0, 0, 0, " ");
+			
+			// timestep
+			uiDefButF(block, NUM, REDRAWBUTSOBJECT, "Timestep:",10,120,100,20, &sphmd->sim_parms->timestep, 0.0, 1.0, 0.0001f, 0, " ");
+			// totaltime
+			uiDefButF(block, NUM, REDRAWBUTSOBJECT, "Totaltime:",110,120,100,20, &sphmd->sim_parms->totaltime, 0.0, 10000.0, 0.001f, 0, " ");
+			uiDefButF(block, NUM, REDRAWBUTSOBJECT, "Sampling distance:",210,120,100,20, &sphmd->sim_parms->samplingdistance, 0.0, 1.0, 0.0001f, 0, " ");
+			uiDefButF(block, NUM, REDRAWBUTSOBJECT, "Smoothing Length:",10,100,100,20, &sphmd->sim_parms->smoothinglength, 0.0, 10.0, 0.1f, 0, " ");
+			
+			uiDefButI(block, NUM, REDRAWBUTSOBJECT, "Surface:",110,100,100,20, &sphmd->sim_parms->computesurfaceevery, 0.0, 1000.0, 1.0, 0, " ");
+			uiDefButI(block, NUM, REDRAWBUTSOBJECT, "Fast Marching:",210,100,100,20, &sphmd->sim_parms->fastmarchingevery, 0.0, 1000.0, 1.0, 0, " ");
+		}
+		else if(sphmd->sim_parms && (sphmd->sim_parms->flags & SPH_SIMSETTINGS_FLAG_FLUID))
+		{
+			uiDefButI(block, NUM, REDRAWBUTSOBJECT, "Resolution:",10,120,100,20, &sphmd->sim_parms->resolution, 1.0, 1000.0, 1.0f, 0, " ");
+		}
+	}
+	
+	
+	uiBlockEndAlign(block);
+}
+
 void object_panels()
 {
 	Object *ob;
@@ -5682,6 +5817,7 @@ void physics_panels()
 		object_panel_cloth(ob);
 		object_panel_cloth_II(ob);
 		object_panel_cloth_III(ob);
+		object_panel_sph(ob);
 		object_panel_fluidsim(ob);
 	}
 }
diff --git a/source/blender/src/drawobject.c b/source/blender/src/drawobject.c
index 561024fb74b..84289be0928 100644
--- a/source/blender/src/drawobject.c
+++ b/source/blender/src/drawobject.c
@@ -1883,11 +1883,13 @@ static void draw_em_measure_stats(Object *ob, EditMesh *em)
 {
 	EditEdge *eed;
 	EditFace *efa;
+	EditVert *eve;
 	float v1[3], v2[3], v3[3], v4[3];
 	float fvec[3];
 	char val[32]; /* Stores the measurement display text here */
 	char conv_float[5]; /* Use a float conversion matching the grid size */
 	float area, col[3]; /* area of the face,  color of the text to draw */
+	int i = 0;
 	
 	/* make the precission of the pronted value proportionate to the gridsize */
 	if ((G.vd->grid) < 0.01)
@@ -1941,7 +1943,7 @@ static void draw_em_measure_stats(Object *ob, EditMesh *em)
 		if(col[1]> 0.5) {col[0]*=0.7; col[2]*= 0.7;}
 		else col[1]= col[1]*0.7 + 0.3;
 		glColor3fv(col);
-		
+		/*
 		for(efa= em->faces.first; efa; efa= efa->next) {
 			if((efa->f & SELECT) || (G.moving && faceselectedOR(efa, SELECT)) ) {
 				VECCOPY(v1, efa->v1->co);
@@ -1966,6 +1968,24 @@ static void draw_em_measure_stats(Object *ob, EditMesh *em)
 				glRasterPos3fv(efa->cent);
 				BMF_DrawString( G.fonts, val);
 			}
+		}*/
+		
+		/* draw IDs of mesh vertexes */
+		for(eve = em->verts.first; eve; eve = eve->next) {
+			char val[32];
+			float fvec[3];
+			VecLerpf(fvec, ob->loc, eve->co, 1.1);
+			glRasterPos3f(fvec[0], fvec[1], fvec[2]);
+
+			sprintf(val, "%d", eve->keyindex);
+			BMF_DrawString(G.fonts, val);
+		}
+		for(efa= em->faces.first; efa; efa= efa->next) {
+			char val[32];
+			sprintf(val, "%d", i);
+			glRasterPos3fv(efa->cent);
+			BMF_DrawString( G.fonts, val);
+			i++;
 		}
 	}