Several minor fixes:

- Added part of Austin's msvc8 fixes (vector::erase function
  was "misused"), hopefully compiles better now.
- Ctrl-b now also bakes a selected fluidsim domain
  similar to the softbodies.
- Added surface smoothing option for domains: default is
  1, higher values result in a smoother surface (and probably
  slightly higher comupation times), while 0 means the surface
  is not modified at all.
- Added BLENDER_ELBEEMBOBJABORT environment variable in readBobj,
  if >0 quits blender when a not yet existing fluidsim
  frame should be loaded. Useful for rendering simulations
  as far as possible from the command line.
- Surface normals pointer is now set to NULL in readfile.c
- Fixed win32 error string handling, now uses a function
  to return the string from the solver.
- Fixed fluidsim particle halo scaling problem.
- Solver update

22 files changed, 518 insertions, 373 deletions

diff --git a/intern/elbeem/extern/LBM_fluidsim.h b/intern/elbeem/extern/LBM_fluidsim.h
index b879dd41396..c93d1292e47 100644
--- a/intern/elbeem/extern/LBM_fluidsim.h
+++ b/intern/elbeem/extern/LBM_fluidsim.h
@@ -68,16 +68,6 @@ int performElbeemSimulation(char *cfgfilename);
 void fluidsimGetAxisAlignedBB(struct Mesh *mesh, float obmat[][4],
 		 /*RET*/ float start[3], /*RET*/ float size[3], /*RET*/ struct Mesh **bbmesh );
 
-// implemented in intern/elbeem/utilities.cpp
-/* set elbeem debug output level (0=off to 10=full on) */
-void elbeemSetDebugLevel(int level);
-/* elbeem debug output function */
-void elbeemDebugOut(char *msg);
-
-/* estimate how much memory a given setup will require */
-double elbeemEstimateMemreq(int res, 
-		float sx, float sy, float sz,
-		int refine, char *retstr);
 
 #endif
 
diff --git a/intern/elbeem/extern/elbeem.h b/intern/elbeem/extern/elbeem.h
index 83a9cb99afc..f0d154c18a4 100644
--- a/intern/elbeem/extern/elbeem.h
+++ b/intern/elbeem/extern/elbeem.h
@@ -73,6 +73,10 @@ typedef struct elbeemSimulationSettings {
 
 	/* global transformation to apply to fluidsim mesh */
 	float surfaceTrafo[4*4];
+
+	/* development variables, testing for upcoming releases...*/
+	float farFieldSize;
+
 } elbeemSimulationSettings;
 
 
@@ -122,8 +126,12 @@ extern "C"  {
 // reset elbeemSimulationSettings struct with defaults
 void elbeemResetSettings(struct elbeemSimulationSettings*);
  
-// start fluidsim init
+// start fluidsim init (returns !=0 upon failure)
 int elbeemInit(struct elbeemSimulationSettings*);
+// get failure message during simulation or init
+// if an error occured (the string is copied into buffer,
+// max. length = 256 chars )
+void elbeemGetErrorString(char *buffer);
 
 // reset elbeemMesh struct with zeroes
 void elbeemResetMesh(struct elbeemMesh*);
@@ -135,17 +143,36 @@ int elbeemAddMesh(struct elbeemMesh*);
 int elbeemSimulate(void);
 
 
-// helper function - simplify animation channels
+// helper functions 
+
+// simplify animation channels
 // returns if the channel and its size changed
 int elbeemSimplifyChannelFloat(float *channel, int *size);
 int elbeemSimplifyChannelVec3(float *channel, int *size);
 
+// helper functions implemented in utilities.cpp
+
+/* set elbeem debug output level (0=off to 10=full on) */
+void elbeemSetDebugLevel(int level);
+/* elbeem debug output function, prints if debug level >0 */
+void elbeemDebugOut(char *msg);
+
+/* estimate how much memory a given setup will require */
+double elbeemEstimateMemreq(int res,
+    float sx, float sy, float sz,
+    int refine, char *retstr);
+
+
+
 #ifdef __cplusplus
 }
 #endif // __cplusplus
 
+
+
 /******************************************************************************/
-// internal defines, do not use for setting up simulation
+// internal defines, do not use for initializing elbeemMesh
+// structs, for these use OB_xxx defines above
 
 /*! fluid geometry init types */
 #define FGI_FLAGSTART   16
diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp
index 8e96879737f..daee5f831d6 100644
--- a/intern/elbeem/intern/elbeem.cpp
+++ b/intern/elbeem/intern/elbeem.cpp
@@ -67,6 +67,8 @@ void elbeemResetSettings(elbeemSimulationSettings *set) {
 	set->generateVertexVectors = 0;
 	set->surfaceSmoothing = 1.;
 
+	set->farFieldSize = 0.;
+
 	// init identity
 	for(int i=0; i<16; i++) set->surfaceTrafo[i] = 0.0;
 	for(int i=0; i<4; i++) set->surfaceTrafo[i*4+i] = 1.0;
@@ -87,6 +89,13 @@ int elbeemInit(elbeemSimulationSettings *settings) {
 	return 0;
 }
 
+// error message access
+extern "C" 
+void elbeemGetErrorString(char *buffer) {
+	if(!buffer) return;
+	strncpy(buffer,gElbeemErrorString,256);
+}
+
 // reset elbeemMesh struct with zeroes
 extern "C" 
 void elbeemResetMesh(elbeemMesh *mesh) {
diff --git a/intern/elbeem/intern/elbeem.h b/intern/elbeem/intern/elbeem.h
index 83a9cb99afc..f0d154c18a4 100644
--- a/intern/elbeem/intern/elbeem.h
+++ b/intern/elbeem/intern/elbeem.h
@@ -73,6 +73,10 @@ typedef struct elbeemSimulationSettings {
 
 	/* global transformation to apply to fluidsim mesh */
 	float surfaceTrafo[4*4];
+
+	/* development variables, testing for upcoming releases...*/
+	float farFieldSize;
+
 } elbeemSimulationSettings;
 
 
@@ -122,8 +126,12 @@ extern "C"  {
 // reset elbeemSimulationSettings struct with defaults
 void elbeemResetSettings(struct elbeemSimulationSettings*);
  
-// start fluidsim init
+// start fluidsim init (returns !=0 upon failure)
 int elbeemInit(struct elbeemSimulationSettings*);
+// get failure message during simulation or init
+// if an error occured (the string is copied into buffer,
+// max. length = 256 chars )
+void elbeemGetErrorString(char *buffer);
 
 // reset elbeemMesh struct with zeroes
 void elbeemResetMesh(struct elbeemMesh*);
@@ -135,17 +143,36 @@ int elbeemAddMesh(struct elbeemMesh*);
 int elbeemSimulate(void);
 
 
-// helper function - simplify animation channels
+// helper functions 
+
+// simplify animation channels
 // returns if the channel and its size changed
 int elbeemSimplifyChannelFloat(float *channel, int *size);
 int elbeemSimplifyChannelVec3(float *channel, int *size);
 
+// helper functions implemented in utilities.cpp
+
+/* set elbeem debug output level (0=off to 10=full on) */
+void elbeemSetDebugLevel(int level);
+/* elbeem debug output function, prints if debug level >0 */
+void elbeemDebugOut(char *msg);
+
+/* estimate how much memory a given setup will require */
+double elbeemEstimateMemreq(int res,
+    float sx, float sy, float sz,
+    int refine, char *retstr);
+
+
+
 #ifdef __cplusplus
 }
 #endif // __cplusplus
 
+
+
 /******************************************************************************/
-// internal defines, do not use for setting up simulation
+// internal defines, do not use for initializing elbeemMesh
+// structs, for these use OB_xxx defines above
 
 /*! fluid geometry init types */
 #define FGI_FLAGSTART   16
diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp
index ed0fab289d8..07bc6b7c855 100644
--- a/intern/elbeem/intern/isosurface.cpp
+++ b/intern/elbeem/intern/isosurface.cpp
@@ -37,7 +37,7 @@ IsoSurface::IsoSurface(double iso) :
   mInitDone(false),
 	mSmoothSurface(0.0), mSmoothNormals(0.0),
 	mAcrossEdge(), mAdjacentFaces(),
-	mCutoff(-1), // off by default
+	mCutoff(-1), mCutArray(NULL),// off by default
 	mFlagCnt(1),
 	mSCrad1(0.), mSCrad2(0.), mSCcenter(0.)
 {
@@ -152,6 +152,7 @@ void IsoSurface::triangulate( void )
 	const int cubieOffsetZ[8] = {
 		0,0,0,0,  1,1,1,1 };
 
+	const int coAdd=2;
   // let the cubes march 
 	pz = mStart[2]-gsz*0.5;
 	for(int k=1;k<(mSizez-2);k++) {
@@ -259,12 +260,17 @@ void IsoSurface::triangulate( void )
 
 				}
 
-				const int coAdd=2;
-				if(i<coAdd+mCutoff) continue;
-				if(j<coAdd+mCutoff) continue;
-				if((mCutoff>0) && (k<coAdd)) continue;
-				if(i>mSizex-2-coAdd-mCutoff) continue;
-				if(j>mSizey-2-coAdd-mCutoff) continue;
+				if( (i<coAdd+mCutoff) ||
+				    (j<coAdd+mCutoff) ||
+				    ((mCutoff>0) && (k<coAdd)) ||// bottom layer
+				    (i>mSizex-2-coAdd-mCutoff) ||
+				    (j>mSizey-2-coAdd-mCutoff) ) {
+					if(mCutArray) {
+						if(k < mCutArray[j*this->mSizex+i]) continue;
+					} else {
+						continue;
+					}
+				}
 
 				// Create the triangles... 
 				for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
diff --git a/intern/elbeem/intern/isosurface.h b/intern/elbeem/intern/isosurface.h
index 90a4e11e95f..5e5115ef236 100644
--- a/intern/elbeem/intern/isosurface.h
+++ b/intern/elbeem/intern/isosurface.h
@@ -95,6 +95,8 @@ class IsoSurface :
 
 		//! cutoff border area
 		int mCutoff;
+		//! cutoff heigh values
+		int *mCutArray;
 		
 		//! trimesh vars
 		vector<int> flags;
@@ -160,6 +162,8 @@ class IsoSurface :
 		}
 		//! set cut off border
 		inline void setCutoff(int set) { mCutoff = set; };
+		//! set cut off border
+		inline void setCutArray(int *set) { mCutArray = set; };
 
 		//! OpenGL viz "interface"
 		unsigned int getIsoVertexCount() {
diff --git a/intern/elbeem/intern/ntl_blenderdumper.cpp b/intern/elbeem/intern/ntl_blenderdumper.cpp
index fce5a085f59..ecd1967c877 100644
--- a/intern/elbeem/intern/ntl_blenderdumper.cpp
+++ b/intern/elbeem/intern/ntl_blenderdumper.cpp
@@ -91,7 +91,7 @@ int ntlBlenderDumper::renderScene( void )
   vector<ntlTriangle> Triangles;
   vector<ntlVec3Gfx>  Vertices;
   vector<ntlVec3Gfx>  VertNormals;
-	errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) );
+	//errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) );
 
 	/* init geometry array, first all standard objects */
 	int idCnt = 0;          // give IDs to objects
@@ -107,13 +107,13 @@ int ntlBlenderDumper::renderScene( void )
 		if(tid & GEOCLASSTID_SHADER) {
 			ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter); //dynamic_cast<ntlGeometryShader*>(*iter);
 			hideObjs.push_back( (*iter)->getName() );
-			geoshad->notifyShaderOfDump(glob->getAniCount(),nrStr,glob->getOutFilename());
+			geoshad->notifyShaderOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename());
 			for (vector<ntlGeometryObject*>::iterator siter = geoshad->getObjectsBegin();
 					siter != geoshad->getObjectsEnd();
 					siter++) {
 				if(debugOut) debMsgStd("ntlBlenderDumper::BuildScene",DM_MSG,"added shader geometry "<<(*siter)->getName(), 8);
 
-				(*siter)->notifyOfDump(glob->getAniCount(),nrStr,glob->getOutFilename());
+				(*siter)->notifyOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename(), this->mSimulationTime);
 				bool doDump = false;
 				bool isPreview = false;
 				// only dump final&preview surface meshes
diff --git a/intern/elbeem/intern/ntl_geometryobject.cpp b/intern/elbeem/intern/ntl_geometryobject.cpp
index 7bbae5d87a1..5ab13a946d3 100644
--- a/intern/elbeem/intern/ntl_geometryobject.cpp
+++ b/intern/elbeem/intern/ntl_geometryobject.cpp
@@ -157,9 +157,9 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
 
 /*! notify object that dump is in progress (e.g. for particles) */
 // default action - do nothing...
-void ntlGeometryObject::notifyOfDump(int frameNr,char *frameNrStr,string outfilename) {
+void ntlGeometryObject::notifyOfDump(int dumtp, int frameNr,char *frameNrStr,string outfilename, double simtime) {
   bool debugOut=false;
-	if(debugOut) debMsgStd("ntlGeometryObject::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<","<<frameNr<<" to "<<outfilename, 10); // DEBUG
+	if(debugOut) debMsgStd("ntlGeometryObject::notifyOfDump",DM_MSG," dt:"<<dumtp<<" obj:"<<this->getName()<<" frame:"<<frameNrStr<<","<<frameNr<<",t"<<simtime<<" to "<<outfilename, 10); // DEBUG
 }
 
 /*****************************************************************************/
diff --git a/intern/elbeem/intern/ntl_geometryobject.h b/intern/elbeem/intern/ntl_geometryobject.h
index 157d160c0e4..ae81d748818 100644
--- a/intern/elbeem/intern/ntl_geometryobject.h
+++ b/intern/elbeem/intern/ntl_geometryobject.h
@@ -16,6 +16,8 @@
 class ntlRenderGlobals;
 class ntlTriangle;
 
+#define DUMP_FULLGEOMETRY 1
+#define DUMP_PARTIAL      2
 
 class ntlGeometryObject : public ntlGeometryClass
 {
@@ -38,7 +40,7 @@ class ntlGeometryObject : public ntlGeometryClass
 				vector<ntlVec3Gfx> *normals, int objectId ) = 0;
 		
 		/*! notify object that dump is in progress (e.g. for particles) */
-		virtual void notifyOfDump(int frameNr,char *frameNrStr,string outfilename);
+		virtual void notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime);
 
 		/*! Search the material for this object from the material list */
 		void searchMaterial(vector<ntlMaterial *> *mat);
diff --git a/intern/elbeem/intern/ntl_geometryshader.h b/intern/elbeem/intern/ntl_geometryshader.h
index 7d767b0f8b9..6adc5629a2d 100644
--- a/intern/elbeem/intern/ntl_geometryshader.h
+++ b/intern/elbeem/intern/ntl_geometryshader.h
@@ -40,7 +40,7 @@ class ntlGeometryShader :
 		virtual vector<ntlGeometryObject *>::iterator getObjectsEnd() { return mObjects.end(); }
 		
 		/*! notify object that dump is in progress (e.g. for field dump) */
-		virtual void notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename) = 0;
+		virtual void notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0;
 
 	protected:
 
diff --git a/intern/elbeem/intern/ntl_world.h b/intern/elbeem/intern/ntl_world.h
index 2256e85ec96..131f1828601 100644
--- a/intern/elbeem/intern/ntl_world.h
+++ b/intern/elbeem/intern/ntl_world.h
@@ -380,9 +380,6 @@ private:
 	bool mSingleFrameMode;
 	//! filename for single frame mode
 	string mSingleFrameFilename;
-
-	/*! Two random number streams for photon generation (one for the directions, the other for russion roulette) */
-	//ntlRandomStream *mpRndDirections, *mpRndRoulette;
 };
 
 
diff --git a/intern/elbeem/intern/particletracer.cpp b/intern/elbeem/intern/particletracer.cpp
index 0a5af90564a..828674675e2 100644
--- a/intern/elbeem/intern/particletracer.cpp
+++ b/intern/elbeem/intern/particletracer.cpp
@@ -20,6 +20,8 @@
 #include <zlib.h>
 
 
+// particle object id counter
+int ParticleObjectIdCnt = 1;
 
 /******************************************************************************
  * Standard constructor
@@ -34,19 +36,21 @@ ParticleTracer::ParticleTracer() :
 	mPartScale(1.0) , mPartHeadDist( 0.5 ), mPartTailDist( -4.5 ), mPartSegments( 4 ),
 	mValueScale(0),
 	mValueCutoffTop(0.0), mValueCutoffBottom(0.0),
-	mDumpParts(0), mShowOnly(0), mpTrafo(NULL)
+	mDumpParts(0), mDumpText(0), mDumpTextFile(""), mShowOnly(0), 
+	mNumInitialParts(0), mpTrafo(NULL)
 {
 	// check env var
-#ifdef ELBEEM_PLUGIN
-	mDumpParts=1; // default on
-#endif // ELBEEM_PLUGIN
-	if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG!
-		int set = atoi( getenv("ELBEEM_DUMPPARTICLE") );
-		if((set>=0)&&(set!=mDumpParts)) {
-			mDumpParts=set;
-			debMsgStd("ParticleTracer::notifyOfDump",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8);
-		}
-	}
+//#ifdef ELBEEM_PLUGIN
+	//mDumpParts=1; // default on
+//#endif // ELBEEM_PLUGIN
+	//if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG!
+		//int set = atoi( getenv("ELBEEM_DUMPPARTICLE") );
+		//if((set>=0)&&(set!=mDumpParts)) {
+			//mDumpParts=set;
+			//debMsgStd("ParticleTrace",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8);
+		//}
+	//}
+	debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10);
 };
 
 ParticleTracer::~ParticleTracer() {
@@ -60,13 +64,9 @@ void ParticleTracer::parseAttrList(AttributeList *att)
 {
 	AttributeList *tempAtt = mpAttrs; 
 	mpAttrs = att;
-	int mNumParticles =0; // UNUSED
-	int mTrailLength  = 0; // UNUSED
-	int mTrailInterval= 0; // UNUSED
-	mNumParticles = mpAttrs->readInt("particles",mNumParticles, "ParticleTracer","mNumParticles", false);
-	mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
-	mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
 
+	mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false);
+	errMsg(" NUMP"," "<<mNumInitialParts);
 	mPartScale    = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false);
 	mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false);
 	mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false);
@@ -75,18 +75,23 @@ void ParticleTracer::parseAttrList(AttributeList *att)
 	mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false);
 	mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false);
 
-	mDumpParts    = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
+	mDumpParts   = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
+	mDumpText    = mpAttrs->readInt  ("part_textdump",mDumpText, "ParticleTracer","mDumpText", false);
 	mShowOnly    = mpAttrs->readInt  ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false);
+	mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false);
 
 	string matPart;
 	matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false);
 	setMaterialName( matPart );
-	// trail length has to be at least one, if anything should be displayed
-	//if((mNumParticles>0)&&(mTrailLength<2)) mTrailLength = 2;
+
+	// unused...
+	int mTrailLength  = 0; // UNUSED
+	int mTrailInterval= 0; // UNUSED
+	mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
+	mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
 
 	// restore old list
 	mpAttrs = tempAtt;
-	//mParts.resize(mTrailLength*mTrailInterval);
 }
 
 /******************************************************************************
@@ -146,7 +151,7 @@ void ParticleTracer::addParticle(float x, float y, float z)
 void ParticleTracer::cleanup() {
 	// cleanup
 	int last = (int)mParts.size()-1;
-	//for(vector<ParticleObject>::iterator pit= getParticlesBegin();pit!= getParticlesEnd(); pit++) {
+	if(mDumpText>0) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }
 
 	for(int i=0; i<=last; i++) {
 		if( mParts[i].getActive()==false ) {
@@ -158,47 +163,18 @@ void ParticleTracer::cleanup() {
 }
 		
 /******************************************************************************
- * save particle positions before adding a new timestep
- * copy "one index up", newest has to remain unmodified, it will be
- * advanced after the next smiulation step
+ *! dump particles if desired 
  *****************************************************************************/
-void ParticleTracer::savePreviousPositions()
-{
-	//debugOut(" PARTS SIZE "<<mParts.size() ,10);
-	/*
-	if(mTrailIntervalCounter==0) {
-	//errMsg("spp"," PARTS SIZE "<<mParts.size() );
-		for(size_t l=mParts.size()-1; l>0; l--) {
-			if( mParts[l].size() != mParts[l-1].size() ) {
-				errFatal("ParticleTracer::savePreviousPositions","Invalid array sizes ["<<l<<"]="<<mParts[l].size()<<
-						" ["<<(l+1)<<"]="<<mParts[l+1].size() <<" , total "<< mParts.size() , SIMWORLD_GENERICERROR);
-				return;
-			}
-
-			for(size_t i=0; i<mParts[l].size(); i++) {
-				mParts[l][i] = mParts[l-1][i];
-			}
-
-		}
-	} 
-	mTrailIntervalCounter++;
-	if(mTrailIntervalCounter>=mTrailInterval) mTrailIntervalCounter = 0;
-	UNUSED!? */
-}
-
-
-/*! dump particles if desired */
-void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilename) {
+void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
 	debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr, 10); // DEBUG
 
-	if(mDumpParts>0) {
+	if(
+			(dumptype==DUMP_FULLGEOMETRY)&&
+			(mDumpParts>0)) {
 		// dump to binary file
 		std::ostringstream boutfilename("");
-		//boutfilename << ecrpath.str() << glob->getOutFilename() <<"_"<< this->getName() <<"_" << frameNrStr << ".obj";
-		//boutfilename << outfilename <<"_particles_"<< this->getName() <<"_" << frameNrStr<< ".gz";
 		boutfilename << outfilename <<"_particles_" << frameNrStr<< ".gz";
-		debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() 
-				<<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
+		debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
 
 		// output to zipped file
 		gzFile gzf;
@@ -210,21 +186,26 @@ void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilenam
 			numParts = 0;
 			for(size_t i=0; i<mParts.size(); i++) {
 				if(!mParts[i].getActive()) continue;
+				//if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
 				numParts++;
 			}
 			gzwrite(gzf, &numParts, sizeof(numParts));
 			for(size_t i=0; i<mParts.size(); i++) {
 				if(!mParts[i].getActive()) { continue; }
-				if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
+				//if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
 				ParticleObject *p = &mParts[i];
-				int type = p->getType();
+				//int type = p->getType();  // export whole type info
+				int type = p->getFlags(); // debug export whole type & status info
 				ntlVec3Gfx pos = p->getPos();
 				float size = p->getSize();
 
 				if(type&PART_FLOAT) { // WARNING same handling for dump!
 					// add one gridcell offset
 					//pos[2] += 1.0; 
-				}
+				} 
+				// display as drop for now externally
+				//else if(type&PART_TRACER) { type |= PART_DROP; }
+
 				pos = (*mpTrafo) * pos;
 
 				ntlVec3Gfx v = p->getVel();
@@ -240,6 +221,64 @@ void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilenam
 			gzclose( gzf );
 		}
 	} // dump?
+
+	// dfor partial & full dump
+	if(mDumpText>0) {
+		// dump to binary file
+		std::ostringstream boutfilename("");
+		if(mDumpTextFile.length()>1) {   
+			boutfilename << mDumpTextFile <<  ".cpart2"; 
+		} else {                           
+			boutfilename << outfilename <<"_particles" <<  ".cpart2"; 
+		}
+		debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
+
+		int numParts = 0;
+		// only dump bubble particles
+		for(size_t i=0; i<mParts.size(); i++) {
+			//if(!mParts[i].getActive()) continue;
+			//if(!(mParts[i].getType()&PART_BUBBLE)) continue;
+			numParts++;
+		}
+
+		// output to text file
+		gzFile gzf;
+		if(frameNr==0) {
+			gzf = gzopen(boutfilename.str().c_str(), "w0");
+
+			gzprintf( gzf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
+			// fixed time scale for now
+			gzprintf( gzf, "T %f \n\n", 1.0);
+		} else {
+			gzf = gzopen(boutfilename.str().c_str(), "a+0");
+		}
+
+		// add current set
+		if(gzf) {
+			gzprintf( gzf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", frameNr, simtime, numParts );
+			gzprintf( gzf, "S %f \n\n", simtime );
+			
+			for(size_t i=0; i<mParts.size(); i++) {
+				ParticleObject *p = &mParts[i];
+				ntlVec3Gfx pos = p->getPos();
+				float size = p->getSize();
+				if(!mParts[i].getActive()) { size=0.; } // switch "off"
+
+				pos = (*mpTrafo) * pos;
+				ntlVec3Gfx v = p->getVel();
+				v[0] *= mpTrafo->value[0][0];
+				v[1] *= mpTrafo->value[1][1];
+				v[2] *= mpTrafo->value[2][2];
+				
+				gzprintf( gzf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
+				gzprintf( gzf, "s %f \n", size );
+				gzprintf( gzf, "\n", size );
+			}
+
+			gzprintf( gzf, "# %d end  ", frameNr );
+			gzclose( gzf );
+		}
+	}
 }
 
 
@@ -293,6 +332,7 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
 			case 2: if(!(type&PART_DROP))   continue; break;
 			case 3: if(!(type&PART_INTER))  continue; break;
 			case 4: if(!(type&PART_FLOAT))  continue; break;
+			case 5: if(!(type&PART_TRACER))  continue; break;
 			}
 		} else {
 			// by default dont display inter
diff --git a/intern/elbeem/intern/particletracer.h b/intern/elbeem/intern/particletracer.h
index b9179b68dc0..453c1a793c0 100644
--- a/intern/elbeem/intern/particletracer.h
+++ b/intern/elbeem/intern/particletracer.h
@@ -16,24 +16,31 @@ template<class Scalar> class ntlMatrix4x4;
 #define PART_DROP   (1<< 2)
 #define PART_INTER  (1<< 3)
 #define PART_FLOAT  (1<< 4)
+#define PART_TRACER (1<< 5)
 
 // particle state
 #define PART_IN     (1<< 8)
 #define PART_OUT    (1<< 9)
 #define PART_INACTIVE (1<<10)
 
+// defines for particle movement
+#define MOVE_FLOATS 1
+#define FLOAT_JITTER 0.03
+
+extern int ParticleObjectIdCnt;
+
 //! A single particle
 class ParticleObject
 {
 	public:
   	//! Standard constructor
   	inline ParticleObject(ntlVec3Gfx mp) :
-			mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0) { };
+			mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0) { mId = ParticleObjectIdCnt++; };
   	//! Copy constructor
   	inline ParticleObject(const ParticleObject &a) :
 			mPos(a.mPos), mVel(a.mVel), mSize(a.mSize), 
 			mStatus(a.mStatus),
-			mLifeTime(a.mLifeTime) { };
+			mLifeTime(a.mLifeTime) { mId = ParticleObjectIdCnt++; };
   	//! Destructor
   	inline ~ParticleObject() { /* empty */ };
 
@@ -81,8 +88,12 @@ class ParticleObject
 		//! set type (lower byte)
 		inline void setLifeTime(int set) { mLifeTime = set; }
 		
+		inline int getId() const { return mId; }
+
 	protected:
 
+		/*! only for debugging */
+		int mId;
 		/*! the particle position */
 		ntlVec3Gfx mPos;
 		/*! the particle velocity */
@@ -109,9 +120,6 @@ class ParticleTracer :
 		//! add a particle at this position
 		void addParticle(float x, float y, float z);
 
-		//! save particle positions before adding a new timestep
-		void savePreviousPositions();
-
 		//! draw the particle array
 		void draw();
 		
@@ -125,6 +133,9 @@ class ParticleTracer :
 		
 		//! get the number of particles
 		inline int  getNumParticles() 				{ return mParts.size(); }
+		//! set/get the number of particles
+		inline void setNumInitialParticles(int set) { mNumInitialParts=set; }
+		inline int  getNumInitialParticles() 	      { return mNumInitialParts; }
 
 		//! iterate over all newest particles (for advancing positions)
 		inline vector<ParticleObject>::iterator getParticlesBegin() { return mParts.begin(); }
@@ -149,7 +160,13 @@ class ParticleTracer :
 		
 		/*! set/get dump flag */
 		inline void setDumpParts(bool set) { mDumpParts = set; }
-		inline bool getDumpParts() { return mDumpParts; }
+		inline bool getDumpParts()         { return mDumpParts; }
+		/*! set/get dump flag */
+		inline void setDumpText(bool set) { mDumpText = set; }
+		inline bool getDumpText()         { return mDumpText; }
+		/*! set/get dump text file */
+		inline void setDumpTextFile(std::string set) { mDumpTextFile = set; }
+		inline std::string getDumpTextFile()         { return mDumpTextFile; }
 		
 		//! set the particle scaling factor
 		inline void setPartScale(float set) { mPartScale = set; }
@@ -161,7 +178,7 @@ class ParticleTracer :
 				vector<ntlVec3Gfx> *vertices, 
 				vector<ntlVec3Gfx> *normals, int objectId );
 
-		virtual void notifyOfDump(int frameNr,char *frameNrStr,string outfilename);
+		virtual void notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename,double simtime);
 		// free deleted particles
 		void cleanup();
 
@@ -194,8 +211,14 @@ class ParticleTracer :
 
 		/*! dump particles (or certain types of) to disk? */
 		int mDumpParts;
+		/*! dump particles (or certain types of) to disk? */
+		int mDumpText;
+		/*! text dump output file */
+		std::string mDumpTextFile;
 		/*! show only a certain type (debugging) */
 		int mShowOnly;
+		/*! no. of particles to init */
+		int mNumInitialParts;
 
 		//! transform matrix
 		ntlMatrix4x4<gfxReal> *mpTrafo;
diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp
index 481eb836ad6..116d68d6824 100644
--- a/intern/elbeem/intern/simulation_object.cpp
+++ b/intern/elbeem/intern/simulation_object.cpp
@@ -156,6 +156,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 	mpLbm->setGeoEnd( mGeoEnd );
 	mpLbm->setRenderGlobals( mpGlob );
 	mpLbm->setName( getName() + "_lbm" );
+	mpLbm->setParticleTracer( mpParts );
 	if(mpElbeemSettings) {
 		// set further settings from API struct init
 		mpLbm->setSmoothing(1.0 * mpElbeemSettings->surfaceSmoothing, 1.0 * mpElbeemSettings->surfaceSmoothing);
@@ -173,6 +174,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 		mpLbm->setDomainBound(dinitType);
 		mpLbm->setDomainPartSlip(mpElbeemSettings->obstaclePartslip);
 		mpLbm->setDumpVelocities(mpElbeemSettings->generateVertexVectors);
+		mpLbm->setFarFieldSize(mpElbeemSettings->farFieldSize);
 		debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<<dinitType<<" ps="<<mpElbeemSettings->obstaclePartslip<<" vv"<<mpElbeemSettings->generateVertexVectors<<","<<mpLbm->getDumpVelocities(), 9 );
 
 		debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<<mpLbm->getGenerateParticles(),10);
@@ -236,7 +238,6 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 	mpParts->setCastShadows( false );
 	mpParts->setReceiveShadows( false );
 	mpParts->searchMaterial( glob->getMaterials() );
-	mpLbm->initParticles(mpParts);
 
 	// this has to be inited here - before, the values might be unknown
 	ntlGeometryObject *surf = mpLbm->getSurfaceGeoObj();
@@ -288,9 +289,6 @@ void SimulationObject::step( void )
 	if(mpParam->getCurrentAniFrameTime()>0.0) {
 		// dont advance for stopped time
 		mpLbm->step();
-
-		mpParts->savePreviousPositions();
-		mpLbm->advanceParticles(mpParts);
 		mTime += mpParam->getTimestep();
 	}
 	if(mpLbm->getPanic()) mPanic = true;
@@ -399,8 +397,8 @@ void SimulationObject::setMouseClick()
 }
 
 /*! notify object that dump is in progress (e.g. for field dump) */
-void SimulationObject::notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename) {
+void SimulationObject::notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
 	if(!mpLbm) return;
-	mpLbm->notifySolverOfDump(frameNr,frameNrStr,outfilename);
+	mpLbm->notifySolverOfDump(dumptype, frameNr,frameNrStr,outfilename);
 }
 
diff --git a/intern/elbeem/intern/simulation_object.h b/intern/elbeem/intern/simulation_object.h
index 4e7bc19294d..b8d1e90596d 100644
--- a/intern/elbeem/intern/simulation_object.h
+++ b/intern/elbeem/intern/simulation_object.h
@@ -91,7 +91,7 @@ class SimulationObject :
 		virtual int postGeoConstrInit(ntlRenderGlobals *glob) { return initializeLbmSimulation(glob); };
 		virtual int initializeShader() { /* ... */ return true; };
 		/*! notify object that dump is in progress (e.g. for field dump) */
-		virtual void notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename);
+		virtual void notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
 		/*! simluation interface: draw the simulation with OpenGL */
 		virtual void draw( void ) {};
 		virtual vector<ntlGeometryObject *>::iterator getObjectsBegin();
diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h
index 7217964e7ba..b24183411ef 100644
--- a/intern/elbeem/intern/solver_class.h
+++ b/intern/elbeem/intern/solver_class.h
@@ -211,7 +211,7 @@ class LbmFsgrSolver :
 		//! finish the init with config file values (allocate arrays...) 
 		virtual bool initializeSolver(); 
 		//! notify object that dump is in progress (e.g. for field dump) 
-		virtual void notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename);
+		virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
 
 #if LBM_USE_GUI==1
 		//! show simulation info (implement LbmSolverInterface pure virtual func)
@@ -267,9 +267,9 @@ class LbmFsgrSolver :
 		/* simulation object interface, just calls stepMain */
 		virtual void step();
 		/*! init particle positions */
-		virtual int initParticles(ParticleTracer *partt);
+		virtual int initParticles();
 		/*! move all particles */
-		virtual void advanceParticles(ParticleTracer *partt );
+		virtual void advanceParticles();
 
 
 		/*! debug object display (used e.g. for preview surface) */
@@ -436,8 +436,6 @@ class LbmFsgrSolver :
 		int mForceTadapRefine;
 		//! border cutoff value
 		int mCutoff;
-		//! store particle tracer
-		ParticleTracer *mpParticles;
 
 		// strict debug interface
 #		if FSGR_STRICT_DEBUG==1
@@ -460,11 +458,10 @@ class LbmFsgrSolver :
 		void initTestdata();
 		void destroyTestdata();
 		void handleTestdata();
-		void exportTestdata();
 		void set3dHeight(int ,int );
 	public:
 		// needed from testdata
-		void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat &retux, LbmFloat &retuy);
+		void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy);
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 
 	public: // former LbmModelLBGK  functions
@@ -700,9 +697,8 @@ class LbmFsgrSolver :
 #define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set) ]
 #define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set) ]
 
-// array data layouts
-// standard array layout  -----------------------------------------------------------------------------------------------
-#define ALSTRING "Standard Array Layout"
+// array handling  -----------------------------------------------------------------------------------------------
+
 //#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
 #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
 #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp
index abe9ef6994a..a417c33caa9 100644
--- a/intern/elbeem/intern/solver_init.cpp
+++ b/intern/elbeem/intern/solver_init.cpp
@@ -396,8 +396,6 @@ LbmFsgrSolver::LbmFsgrSolver() :
 		mGaussw[n] = mGaussw[n]/totGaussw;
 	}
 
-	mpParticles = NULL;
-	//addDrop(false,0,0);
 }
 
 /*****************************************************************************/
@@ -486,6 +484,8 @@ void LbmFsgrSolver::parseAttrList()
 #else // LBM_INCLUDE_TESTSOLVERS!=1
 	// off by default
 	mUseTestdata = 0;
+	if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
+	if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
 #endif // LBM_INCLUDE_TESTSOLVERS!=1
 }
 
@@ -573,7 +573,14 @@ void LbmFsgrSolver::initLevelOmegas()
  *****************************************************************************/
 bool LbmFsgrSolver::initializeSolver()
 {
-//  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... (Layout:"<<ALSTRING<<") "<<this->mInitDone<<" "<<((int)this),1);
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<this->mInitDone<<" "<<(void*)this,1);
+
+	// init cppf stage
+	if(mCppfStage>0) {
+		this->mSizex *= mCppfStage;
+		this->mSizey *= mCppfStage;
+		this->mSizez *= mCppfStage;
+	}
 
 	// size inits to force cubic cells and mult4 level dimensions
 	// and make sure we dont allocate too much...
@@ -584,9 +591,6 @@ bool LbmFsgrSolver::initializeSolver()
 	double sizeReduction = 1.0;
 	double memEstFromFunc = -1.0;
 	string memreqStr("");	
-#if LBM_INCLUDE_TESTSOLVERS==1
-	if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
-#endif
 	while(!memOk) {
 		initGridSizes( this->mSizex, this->mSizey, this->mSizez,
 				this->mvGeoStart, this->mvGeoEnd, mMaxRefine, PARALLEL);
@@ -750,16 +754,25 @@ bool LbmFsgrSolver::initializeSolver()
 	mMaxTimestep = this->mpParam->getTimestep();
 
 	// init isosurf
+	this->mpIso->setIsolevel( this->mIsoValue );
 #if LBM_INCLUDE_TESTSOLVERS==1
 	if(mUseTestdata) {
 		mpTest->setMaterialName( this->mpIso->getMaterialName() );
 		delete this->mpIso;
 		this->mpIso = mpTest;
+		if(mpTest->mDebugvalue1>0.0) { // 3d off
+			mpTest->setIsolevel(-100.0);
+		} else {
+			mpTest->setIsolevel( this->mIsoValue );
+		}
 	}
 #endif // ELBEEM_PLUGIN!=1
-	this->mpIso->setIsolevel( this->mIsoValue );
 	// approximate feature size with mesh resolution
 	float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
+	// smooth vars defined in solver_interface, set by simulation object
+	// reset for invalid values...
+	if((this->mSmoothSurface<0.)||(this->mSmoothSurface>50.)) this->mSmoothSurface = 1.;
+	if((this->mSmoothNormals<0.)||(this->mSmoothNormals>50.)) this->mSmoothNormals = 1.;
 	this->mpIso->setSmoothSurface( this->mSmoothSurface * featureSize );
 	this->mpIso->setSmoothNormals( this->mSmoothNormals * featureSize );
 
@@ -984,6 +997,12 @@ bool LbmFsgrSolver::initializeSolver()
 			errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
 			this->mOutputSurfacePreview = 0; }
 	}
+#if LBM_USE_GUI==1
+	if(this->mOutputSurfacePreview) {
+		errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
+		this->mOutputSurfacePreview = 0; }
+#endif // LBM_USE_GUI==1
+	
 	if(this->mOutputSurfacePreview) {
 
 		// same as normal one, but use reduced size
@@ -1017,13 +1036,17 @@ bool LbmFsgrSolver::initializeSolver()
 
 
 	// now really done...
-  debugOut("LbmFsgrSolver::initialize : Init done ...",10);
+  debugOut("LbmFsgrSolver::initialize : SurfaceGen: SmsOrg("<<this->mSmoothSurface<<","<<this->mSmoothNormals<<","<<featureSize<<"), Iso("<<this->mpIso->getSmoothSurface()<<","<<this->mpIso->getSmoothNormals()<<") ",10);
+  debugOut("LbmFsgrSolver::initialize : Init done ... ",10);
 	this->mInitDone = 1;
 
 #if LBM_INCLUDE_TESTSOLVERS==1
 	initTestdata();
 #endif // ELBEEM_PLUGIN!=1
+	// not inited? dont use...
+	if(mCutoff<0) mCutoff=0;
 
+	initParticles();
 	return true;
 }
 
diff --git a/intern/elbeem/intern/solver_interface.cpp b/intern/elbeem/intern/solver_interface.cpp
index 04e0456286b..ee2f8f281be 100644
--- a/intern/elbeem/intern/solver_interface.cpp
+++ b/intern/elbeem/intern/solver_interface.cpp
@@ -34,7 +34,7 @@ LbmSolverInterface::LbmSolverInterface() :
 	mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop(  (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
   mInitDone( false ),
   mInitDensityGradient( false ),
-	mpAttrs( NULL ), mpParam( NULL ),
+	mpAttrs( NULL ), mpParam( NULL ), mpParticles(NULL),
 	mNumParticlesLost(0), 
 	mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
 	mDebugVelScale( 0.01 ), mNodeInfoString("+"),
@@ -53,7 +53,7 @@ LbmSolverInterface::LbmSolverInterface() :
 	mDumpVelocities(false),
 	mMarkedCells(), mMarkedCellIndex(0),
 	mDomainBound("noslip"), mDomainPartSlipValue(0.1),
-	mTForceStrength(0.0)
+	mTForceStrength(0.0), mFarFieldSize(0.), mCppfStage(0)
 {
 #if ELBEEM_PLUGIN==1
 	if(gDebugLevel<=1) mSilent = true;
@@ -242,7 +242,12 @@ void LbmSolverInterface::parseStdAttrList() {
 	
 	// new test vars
 	mTForceStrength = mpAttrs->readFloat("tforcestrength", mTForceStrength,"LbmSolverInterface", "mTForceStrength", false);
+	mFarFieldSize = mpAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false);
+	// old compat
+	float sizeScale = mpAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false);
+	if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); }
 
+	mCppfStage = mpAttrs->readFloat("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false);
 	mPartGenProb = mpAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false);
 }
 
diff --git a/intern/elbeem/intern/solver_interface.h b/intern/elbeem/intern/solver_interface.h
index c86bbce2a89..ee4ecd79077 100644
--- a/intern/elbeem/intern/solver_interface.h
+++ b/intern/elbeem/intern/solver_interface.h
@@ -260,7 +260,7 @@ class LbmSolverInterface
 		virtual bool initializeSolver() =0;
 		
 		/*! notify object that dump is in progress (e.g. for field dump) */
-		virtual void notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename) = 0;
+		virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0;
 
 		/*! parse a boundary flag string */
 		CellFlagType readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed);
@@ -273,8 +273,8 @@ class LbmSolverInterface
 		virtual void prepareVisualization() { /* by default off */ };
 
 		/*! particle handling */
-		virtual int initParticles(ParticleTracer *partt) = 0;
-		virtual void advanceParticles(ParticleTracer *partt ) = 0;
+		virtual int initParticles() = 0;
+		virtual void advanceParticles() = 0;
 		/*! get surface object (NULL if no surface) */
 		ntlGeometryObject* getSurfaceGeoObj() { return mpIso; }
 
@@ -329,6 +329,9 @@ class LbmSolverInterface
 		inline void setParametrizer(Parametrizer *set) { mpParam = set; }
 		/*! get parametrizer pointer */
 		inline Parametrizer *getParametrizer() { return mpParam; }
+		/*! get/set particle pointer */
+		inline void setParticleTracer(ParticleTracer *set) { mpParticles = set; }
+		inline ParticleTracer *getParticleTracer() { return mpParticles; }
 
 		/*! set density gradient init from e.g. init test cases */
 		inline void setInitDensityGradient(bool set) { mInitDensityGradient = set; }
@@ -379,6 +382,12 @@ class LbmSolverInterface
 		//! set/get dump velocities flag
 		inline void setDomainPartSlip(LbmFloat set)	{ mDomainPartSlipValue = set; }
 		inline LbmFloat getDomainPartSlip() const	{ return mDomainPartSlipValue; }
+		//! set/get far field size
+		inline void setFarFieldSize(LbmFloat set)	{ mFarFieldSize = set; }
+		inline LbmFloat getFarFieldSize() const	{ return mFarFieldSize; }
+		//! set/get cp stage
+		inline void setCpStage(int set)	{ mCppfStage = set; }
+		inline int getCpStage() const	{ return mCppfStage; }
 
 		// cell iterator interface
 		
@@ -474,6 +483,8 @@ class LbmSolverInterface
 
 		/*! get parameters from this parametrize in finishInit */
 		Parametrizer *mpParam;
+		//! store particle tracer
+		ParticleTracer *mpParticles;
 
 		/*! number of particles lost so far */
 		int mNumParticlesLost;
@@ -553,6 +564,10 @@ class LbmSolverInterface
 		//! test vars
 		// strength of applied force
 		LbmFloat mTForceStrength;
+		// 
+		LbmFloat mFarFieldSize;
+		// 
+		int mCppfStage;
 
 };
 
diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp
index c8f0678a6f0..10fca60c546 100644
--- a/intern/elbeem/intern/solver_main.cpp
+++ b/intern/elbeem/intern/solver_main.cpp
@@ -108,7 +108,7 @@ void LbmFsgrSolver::stepMain()
 		int avgcls = (int)(mAvgNumUsedCells/(LONGINT)this->mStepCnt);
   	debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime<<
 		//debMsgDirect( 
-			"mlsups(curr:"<<this->mMLSUPS<<
+			" mlsups(curr:"<<this->mMLSUPS<<
 			" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
 			" totcls:"<<(this->mNumUsedCells)<< sepStr<<
 			" avgcls:"<< avgcls<< sepStr<<
@@ -122,9 +122,6 @@ void LbmFsgrSolver::stepMain()
 			" probs:"<<mNumProblems<< sepStr<<
 			" simt:"<<mSimulationTime<< sepStr<<
 			" for '"<<this->mName<<"' " , 10);
-
-		debMsgDirect(std::endl);
-		debMsgDirect(this->mStepCnt<<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<") ");
 		debMsgDirect(std::endl);
 
 		// nicer output
@@ -208,6 +205,9 @@ void LbmFsgrSolver::stepMain()
 	if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); }
 #endif
 
+	// advance positions with current grid
+	advanceParticles();
+
 	// one of the last things to do - adapt timestep
 	// was in fineAdvance before... 
 	if(mTimeAdap) {
@@ -220,6 +220,13 @@ void LbmFsgrSolver::stepMain()
 	if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
 	if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
 #endif // WIN32
+
+	// output total step time
+	timeend = getTime();
+	debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt
+			<<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
+			<<" totst:"<<getTimeString(timeend-timestart), 3);
+
  //#endif // ELBEEM_PLUGIN!=1
 }
 
@@ -283,11 +290,12 @@ LbmFsgrSolver::mainLoop(int lev)
 	int      calcCellsFilled = this->mNumFilledCells;
 	int      calcCellsEmptied = this->mNumEmptiedCells;
 	int      calcNumUsedCells = this->mNumUsedCells;
+	const int cutMin  = 1;
+	const int cutConst = mCutoff+1;
 
 #	if LBM_INCLUDE_TESTSOLVERS==1
-	if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { 
-		// 3d region off... quit
-		this->mpIso->setIsolevel(-100.0); return; }
+	// 3d region off... quit
+	if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; }
 #endif // ELBEEM_PLUGIN!=1
 	//printLbmCell(lev, 6,6,16, mLevel[lev].setCurr ); // DEBUG
 	
@@ -814,13 +822,18 @@ LbmFsgrSolver::mainLoop(int lev)
 		// for inflow no pargen test
 		// NOBUBBB!
 		if((this->mInitDone) //&&(mUseTestdata) 
-				&& (!((oldFlag|newFlag)&CFNoNbEmpty)) 
+				//&& (!((oldFlag|newFlag)&CFNoNbEmpty)) 
 				&& (!((oldFlag|newFlag)&CFNoDelete)) 
 				&& (this->mPartGenProb>0.0)) {
+			bool doAdd = true;
 			bool bndOk=true;
-			if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)||
-					(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)||
-					(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; }
+			//if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)||
+					//(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)||
+					//(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; }
+			if( (i<cutMin)||(i>this->mSizex-cutMin)||
+					(j<cutMin)||(j>this->mSizey-cutMin)||
+					(k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
+			if(!bndOk) doAdd=false;
 			
 			LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
 			LbmFloat rux = (ux * realWorldFac);
@@ -828,12 +841,19 @@ LbmFsgrSolver::mainLoop(int lev)
 			LbmFloat ruz = (uz * realWorldFac);
 			LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
 			// WHMOD
-			const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor!
-			bool doAdd = true;
 
 			LbmFloat prob = (rand()/(RAND_MAX+1.0));
 			LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl;
-			if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>2.5) ) {
+
+			// reduce probability in outer region
+			if( (i<cutConst)||(i>this->mSizex-cutConst) ){ prob *= 0.5; }
+			if( (j<cutConst)||(j>this->mSizey-cutConst) ){ prob *= 0.5; }
+			if( (k<cutConst)||(k>this->mSizez-cutConst) ){ prob *= 0.5; }
+
+//#define RWVEL_THRESH 1.0
+#define RWVEL_THRESH 1.5
+#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
+			if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
 				// add
 			} else {
 				doAdd = false; // dont...
@@ -851,12 +871,12 @@ LbmFsgrSolver::mainLoop(int lev)
 
 			// "wind" disturbance
 			// use realworld relative velocity here instead?
-			if( 
-					((rl>1.0) && (lcsmqo<P_LCSMQO)) // normal checks
-					||(k>this->mSizez-SLOWDOWNREGION) ) {
+			if( (doAdd) && (
+					((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) // normal checks
+					||(k>this->mSizez-SLOWDOWNREGION)   )) {
 				LbmFloat nuz = uz;
 				LbmFloat jdf; // = 0.05 * (rand()/(RAND_MAX+1.0));
-				if(rl>1.0) jdf *= (rl-1.0);
+				if(rl>RWVEL_WINDTHRESH) jdf *= (rl-RWVEL_WINDTHRESH);
 				if(k>this->mSizez-SLOWDOWNREGION) {
 					// special case
 					LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) );
@@ -881,7 +901,6 @@ LbmFsgrSolver::mainLoop(int lev)
 			if(usqr<0.0001) doAdd=false;   // TODO check!?
 			// if outside, and 20% above sea level, delete, TODO really check level?
 			//if((!bndOk)&&((LbmFloat)k>pTest->mFluidHeight*1.5)) { doAdd=true; } // FORCEDISSOLVE
-			//? if(!bndOk) doAdd=false;
 			//if(this->mStepCnt>700) errMsg("DFJITT"," at "<<PRINT_IJK<<"rwl:"<<rl<<"  usqr:"<<usqr <<" qo:"<<lcsmqo<<" add="<<doAdd );
 
 			if( (doAdd)  ) { // ADD DROP
@@ -911,15 +930,14 @@ LbmFsgrSolver::mainLoop(int lev)
 				LbmFloat srcj = j+0.5+jpy; // TEST? + (pv[1]*1.41);
 				LbmFloat srck = k+0.5+jpz; // TEST? + (pv[2]*1.41);
 				int type=0;
-				//if((s%3)!=2) {
-
-					type=PART_DROP;
-					// drop
-					srci += (pv[0]*1.41);
-					srcj += (pv[1]*1.41);
-					srck += (pv[2]*1.41);
-					if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction
-				//} else { type=PART_FLOAT; }
+				//if((s%3)!=2) {} else { type=PART_FLOAT; }
+				//type = PART_DROP;
+				type = PART_INTER;
+				// drop
+				/*srci += (pv[0]*1.41);
+				srcj += (pv[1]*1.41);
+				srck += (pv[2]*1.41);
+				if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction */
 
 				pv *= len;
 				LbmFloat size = 1.0+ 9.0* (rand()/(RAND_MAX+1.0));
@@ -929,18 +947,20 @@ LbmFsgrSolver::mainLoop(int lev)
 				//? mpParticles->getLast()->advanceVel(); // advance a bit outwards
 				mpParticles->getLast()->setStatus(PART_IN);
 				mpParticles->getLast()->setType(type);
-				//mpParticles->getLast()->setType(PART_INTER);
 				//if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
 				mpParticles->getLast()->setSize(size);
 				//errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<PRINT_VEC(ux,uy,uz) <<" RWu="<<PRINT_VEC(rux,ruy,ruz)<<" add"<<doAdd<<" pvel="<<pv );
 				//errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
 				//if(!bndOk){ mass -= val2fac*size*0.02; } // FORCEDISSOLVE
-				mass -= val2fac*size*0.0015; // NTEST!
+				//const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor!
+				//mass -= val2fac*size*0.0015; // NTEST!
 				//mass -= val2fac*size*0.001; // NTEST!
+				mass -= size*0.0020; // NTEST!
 #if LBMDIM==2
 				mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
 				mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
 #endif // LBMDIM==2
+    		//errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel()  );
 				} // multiple parts
 			} // doAdd
 		} // */
@@ -2487,17 +2507,19 @@ void LbmFsgrSolver::reinitFlags( int workSet ) {
 
 	/* remove empty interface cells that are not allowed to be removed anyway
 	 * this is important, otherwise the dreaded cell-type-flickering can occur! */
-  for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
-       iter != mListEmpty.end(); iter++ ) {
-    int i=iter->x, j=iter->y, k=iter->z;
+  //for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); iter != mListEmpty.end(); iter++ ) {
+  //int i=iter->x, j=iter->y, k=iter->z;
+	//iter = mListEmpty.erase(iter); iter--; // and continue with next...
+  for(int n=0; n<(int)mListEmpty.size(); n++) {
+    int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
 		if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
-			// remove entry
-			if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
-			iter = mListEmpty.erase(iter); 
-			iter--; // and continue with next...
-
 			// treat as "new inter"
 			addToNewInterList(i,j,k);
+			// remove entry
+			if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
+			if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
+			mListEmpty.pop_back();
+			n--; 
 		}
 	} 
 
diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h
index 589f2197422..605bc02c0cf 100644
--- a/intern/elbeem/intern/solver_relax.h
+++ b/intern/elbeem/intern/solver_relax.h
@@ -17,102 +17,9 @@
 
 
 
-#if LBM_INCLUDE_TESTSOLVERS!=1
-
 // off for non testing
 #define PRECOLLIDE_MODS(rho,ux,uy,uz) 
 
-#else // LBM_INCLUDE_TESTSOLVERS!=1
-
-#define PRECOLLIDE_GETPOS \
-	LbmVec( \
-			((this->mvGeoEnd[0]-this->mvGeoStart[0])/(LbmFloat)this->mSizex) * (LbmFloat)i + this->mvGeoStart[0], \
-			((this->mvGeoEnd[1]-this->mvGeoStart[1])/(LbmFloat)this->mSizey) * (LbmFloat)j + this->mvGeoStart[1], \
-			((this->mvGeoEnd[2]-this->mvGeoStart[2])/(LbmFloat)this->mSizez) * (LbmFloat)k + this->mvGeoStart[2]  \
-			)
-
-// debug modifications of collide vars (testing)
-#define PRECOLLIDE_MODS(rho,ux,uy,uz) \
-	if( (this->mTForceStrength>0.) && (RFLAG(0,i,j,k, mLevel[0].setCurr)&CFNoBndFluid) ) { \
-		LbmVec pos = PRECOLLIDE_GETPOS; \
-		LbmVec vel(ux,uy,uz); \
-		mpTest->mControlParts.modifyVelocity(pos,vel); /* */\
-		if((i==16)&&(j==10)){ /*debugMarkCell(0,16,10,0);*/ errMsg("FTDEB"," at "<<PRINT_IJK<<" targ:"<<vel<<",len:"<<norm(vel)<<",  org:"<<ntlVec3Gfx(ux,uy,uz) ); }\
-		ux = vel[0]; uy=vel[1]; uz=vel[2]; \
-		/* test acutal values...? */ \
-	}
- 	/* end PRECOLLIDE_MODS */
-
-// debug modifications of collide vars (testing)
-#define _PRECOLLIDE_MODS(rho,ux,uy,uz) \
-	if(this->mTForceStrength>0.) { \
-		LbmVec u(ux,uy,uz); \
-		LbmVec pos = PRECOLLIDE_GETPOS; \
-		LbmVec targv = u; const int lev=0; \
-		mpTest->mControlParts.modifyVelocity(pos,targv); /* */\
-		LbmVec devia = (targv-u); \
-		LbmVec set; \
-		set = u + (devia/this->mOmega) * this->mTForceStrength; /* */\
-		set = u + targv*this->mTForceStrength; /* */\
-		ux = set[0]; uy=set[1]; uz=set[2]; \
-		if(j==10){ errMsg("FTDEB"," at "<<PRINT_IJK<<" targ"<<targv<<","<<norm(targv)<<"  set:"<<set<<","<<norm(set) ); }\
-	}
- 	/* end PRECOLLIDE_MODS */
-
-#endif // LBM_INCLUDE_TESTSOLVERS!=1
-
-/*
-
-	\
-	if((0) && (j>=0.01*this->mSizey) && (j<0.9*this->mSizey)) { \
-		if((1) && (i>=0.5*this->mSizex) && (i<0.6*this->mSizex)) { \
-			LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \
-			LbmVec u(ux,uy,uz); LbmVec t(1., 0., 0.); \
-			ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; \
-		} \
-		if((1) && (i>=0.65*this->mSizex) && (i<0.75*this->mSizex)) { \
-			LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \
-			LbmVec u(ux,uy,uz); LbmVec t(0., 1., 0.); \
-			ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; \
-		} \
-	} \
-	\
-	if((0) && (j>=0.0*this->mSizey) && (j<1.0 *this->mSizey)) { \
-		if((1) && (i>=0.6 *this->mSizex) && (i<0.7 *this->mSizex)) { \
-			LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \
-			LbmVec u(ux,uy,uz); \
-			LbmVec t(0., 1., 0.); \
-			LbmFloat devia = norm(u-t); \
-			// 0.0001-strong, 0.000001-small, 
-			t = u - (devia/this->mOmega) * this->mTForceStrength; \
-			// ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; 
-			ux = t[0]; uy=t[1]; uz=t[2]; \
-			//ux= uy= uz= 0.0; 
-			//errMsg("DDDD"," at "<<PRINT_IJK);  
-		} \
-	} \
-	if(0) { \
-		LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \
-		LbmVec u(ux,uy,uz); \
-		LbmVec p(i/(LbmFloat)this->mSizex, j/(LbmFloat)this->mSizey, k/(LbmFloat)this->mSizez); \
-		LbmVec cp(0.5, 0.5, 0.5); \
-		LbmVec delt = cp - p; \
-		LbmFloat dist = norm(delt); \
-		normalize(delt); \
-		LbmVec tang = cross(delt, LbmVec(0.,0.,1.)); \
-		normalize(tang); \
-		const LbmFloat falloff = 5.0; \
-		LbmVec targv = tang*1.0 * 1.0/(1.0+falloff*dist); \
-		LbmVec devia = (targv-u); \
-		LbmFloat devial = norm(devia); \
-		LbmVec set; \
-		set = u +targv * this->mTForceStrength; \
-		set = u + (devia/this->mOmega) * this->mTForceStrength; \
-		ux = set[0]; uy=set[1]; uz=set[2]; \
-		if(j==10){ errMsg("FTDEB"," at "<<PRINT_IJK<<" dist:"<<delt<<","<<dist<<" targ"<<targv<<","<<norm(targv)<<"  set:"<<set<<","<<norm(set) ); }\
-	} \
-
-*/
 	
 /******************************************************************************
  * normal relaxation
@@ -1152,14 +1059,14 @@ inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmF
 }
 
 #define DEBUG_CALCPRINTCELL(str,df) {\
-		LbmFloat rho=df[0], ux=0., uy=0., uz=0.; \
+		LbmFloat prho=df[0], pux=0., puy=0., puz=0.; \
 		for(int l=1; l<this->cDfNum; l++) { \
-			rho += df[l];  \
-			ux  += (this->dfDvecX[l]*df[l]);  \
-			uy  += (this->dfDvecY[l]*df[l]);  \
-			uz  += (this->dfDvecZ[l]*df[l]);  \
+			prho += df[l];  \
+			pux  += (this->dfDvecX[l]*df[l]);  \
+			puy  += (this->dfDvecY[l]*df[l]);  \
+			puz  += (this->dfDvecZ[l]*df[l]);  \
 		} \
-		errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<rho<<" vel="<<ntlVec3Gfx(ux,uy,uz) ); \
+		errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<prho<<" vel="<<ntlVec3Gfx(pux,puy,puz) ); \
 	} /* END DEBUG_CALCPRINTCELL */ 
 
 // "normal" collision
diff --git a/intern/elbeem/intern/solver_util.cpp b/intern/elbeem/intern/solver_util.cpp
index 7406c55acbc..8cd698eb65f 100644
--- a/intern/elbeem/intern/solver_util.cpp
+++ b/intern/elbeem/intern/solver_util.cpp
@@ -110,39 +110,6 @@ void LbmFsgrSolver::prepareVisualization( void ) {
 		*this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
 	}
 
-	
-#if LBM_INCLUDE_TESTSOLVERS==1
-	/*if(mUseTestdata) {
-		int border = 1;
-		for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++)
-			for(int j=0;j<mLevel[mMaxRefine].lSizey-1;j++) {
-				for(int l=0; l<=border; l++) {
-					*this->mpIso->lbmGetData( l-1,                         j,ZKOFF) = *this->mpIso->lbmGetData( border+1,                         j,ZKOFF);
-					*this->mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-l, j,ZKOFF) = *this->mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-border-1, j,ZKOFF);
-				}
-			}
-
-		for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++)
-			for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) {      
-				for(int l=0; l<=border; l++) {
-					*this->mpIso->lbmGetData( i, l-1,                         ZKOFF) = *this->mpIso->lbmGetData( i, border+1,                         ZKOFF);
-					*this->mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-l, ZKOFF) = *this->mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-border-1, ZKOFF);
-				}
-			}
-
-		if(LBMDIM == 3) {
-			// only for 3D
-			for(int j=-1;j<mLevel[mMaxRefine].lSizey+1;j++)
-				for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) {      
-					for(int l=0; l<=border; l++) {
-						*this->mpIso->lbmGetData( i,j,l-1                        ) = *this->mpIso->lbmGetData( i,j, border+1                         );
-						*this->mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-l) = *this->mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-1-border);
-					}
-				}
-		}
-	} // testdata */
-#endif // LBM_INCLUDE_TESTSOLVERS==1
-	// */
 
 	// update preview, remove 2d?
 	if((this->mOutputSurfacePreview)&&(LBMDIM==3)) {
@@ -175,7 +142,7 @@ void LbmFsgrSolver::prepareVisualization( void ) {
 				*mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , this->mSizez-1);
 			}
 
-		if(mUseTestdata) {
+		if(mFarFieldSize>=2.) {
 			// also remove preview border
 			for(int k= 0; k< (pvsz-1); ++k) {
 				for(int j=0;j< pvsy;j++) {
@@ -267,11 +234,11 @@ vector<ntlGeometryObject*> LbmFsgrSolver::getDebugObjects() {
  *****************************************************************************/
 
 /*! init particle positions */
-int LbmFsgrSolver::initParticles(ParticleTracer *partt) { 
+int LbmFsgrSolver::initParticles() { 
   int workSet = mLevel[mMaxRefine].setCurr;
   int tries = 0;
   int num = 0;
-	mpParticles=partt;
+	ParticleTracer *partt = mpParticles;
 
   partt->setStart( this->mvGeoStart + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
   partt->setEnd  ( this->mvGeoEnd   + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
@@ -279,11 +246,14 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) {
   partt->setSimStart( ntlVec3Gfx(0.0) );
   partt->setSimEnd  ( ntlVec3Gfx(this->mSizex,   this->mSizey,   getForZMaxBnd(mMaxRefine)) );
   
-  while( (num<partt->getNumParticles()) && (tries<100*partt->getNumParticles()) ) {
+  while( (num<partt->getNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) {
     LbmFloat x,y,z;
-    x = 0.0+(( (LbmFloat)(this->mSizex-1) )     * (rand()/(RAND_MAX+1.0)) );
-    y = 0.0+(( (LbmFloat)(this->mSizey-1) )     * (rand()/(RAND_MAX+1.0)) );
-    z = 0.0+(( (LbmFloat) getForZMax1(mMaxRefine) )* (rand()/(RAND_MAX+1.0)) );
+    //x = 0.0+(( (LbmFloat)(this->mSizex-1) )     * (rand()/(RAND_MAX+1.0)) );
+    //y = 0.0+(( (LbmFloat)(this->mSizey-1) )     * (rand()/(RAND_MAX+1.0)) );
+    //z = 0.0+(( (LbmFloat) getForZMax1(mMaxRefine) )* (rand()/(RAND_MAX+1.0)) );
+    x = 1.0+(( (LbmFloat)(this->mSizex-3.) )          * (rand()/(RAND_MAX+1.0)) );
+    y = 1.0+(( (LbmFloat)(this->mSizey-3.) )          * (rand()/(RAND_MAX+1.0)) );
+    z = 1.0+(( (LbmFloat) getForZMax1(mMaxRefine)-2. )* (rand()/(RAND_MAX+1.0)) );
     int i = (int)(x+0.5);
     int j = (int)(y+0.5);
     int k = (int)(z+0.5);
@@ -292,12 +262,19 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) {
       z = 0.5; // place in the middle of domain
     }
 
-    if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ||
-        TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { // only fluid cells?
+    //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ||
+        //TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { // only fluid cells?
+    if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) 
+        //&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid ) 
+				) { // inner fluid only 
+			bool cellOk = true;
+			FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; }
+			if(!cellOk) continue;
       // in fluid...
       partt->addParticle(x,y,z);
 			partt->getLast()->setStatus(PART_IN);
-			partt->getLast()->setType(PART_BUBBLE);
+			partt->getLast()->setType(PART_TRACER);
+			partt->getLast()->setSize(1.0);
       num++;
     }
     tries++;
@@ -337,6 +314,7 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) {
 				if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
 			}
 		}
+		// place floats on rectangular region FLOAT_JITTER_BND
 		if(mpTest->mDebugvalue2==-10.0){ 
 			const int lev = mMaxRefine;
 			const int sx = mLevel[lev].lSizex;
@@ -348,8 +326,9 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) {
 					LbmFloat x,y,z;
 					x = 0.0+(LbmFloat)(i);
 					y = 0.0+(LbmFloat)(j);
-					//z = 0.5+(LbmFloat)(k);
-					z = mLevel[lev].lSizez/20.0*8.0 - 1.0;
+					//z = mLevel[lev].lSizez/10.0*2.5 - 1.0;
+					z = mLevel[lev].lSizez/20.0*7.5 - 1.0;
+					//z = mLevel[lev].lSizez/20.0*4.5 - 1.0;
 					partt->addParticle(x,y,z);
 					//if( (i>0)&&(i<sx) && (j>0)&&(j<sy) ) { partt->getLast()->setStatus(PART_IN); } else { partt->getLast()->setStatus(PART_OUT); }
 					partt->getLast()->setStatus(PART_IN);
@@ -363,17 +342,21 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) {
 #endif // LBM_INCLUDE_TESTSOLVERS
 
 	
-  debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb, 10);
+  debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb<<", tries:"<<tries, 10);
   if(num != partt->getNumParticles()) return 1;
 
 	return 0;
 }
 
-void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { 
+#define P_CHANGETYPE(p, newtype) \
+		p->setLifeTime(0); \
+    /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \
+		p->setType(newtype); 
+
+void LbmFsgrSolver::advanceParticles() { 
   int workSet = mLevel[mMaxRefine].setCurr;
 	LbmFloat vx=0.0,vy=0.0,vz=0.0;
 	LbmFloat rho, df[27]; //feq[27];
-	if(mpParticles!=partt) { errMsg("LbmFsgrSolver::advanceParticles","Invalid ParticleTracer..."); }
 #define DEL_PART { \
 	/*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<"  ");  */ \
 	p->setActive( false ); \
@@ -394,6 +377,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 	const LbmFloat v1 = 9.0; // v max
 	const LbmFloat v2 = 2.0; // v min
 	const LbmVec rwgrav = vec2L( this->mpParam->getGravity(mSimulationTime) );
+	const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){ 
 
 	// TODO use timestep size
 	//bool isIn,isOut,isInZ;
@@ -402,11 +386,11 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 	//const int cutval = 0; // TODO FIXME add half border!
 	const int cutval = mCutoff; // use full border!?
 	int actCnt=0;
-	if(this->mStepCnt%50==49) { partt->cleanup(); }
-  for(vector<ParticleObject>::iterator pit= partt->getParticlesBegin();
-      pit!= partt->getParticlesEnd(); pit++) {
-    //errMsg("PIT"," pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel()  );
-    //errMsg("PIT"," pit pos:"<< (*pit).getPos()<<" vel:"<< (*pit).getVel()<<" status:"<<convertFlags2String((*pit).getFlags()) <<" " <<partt->getStart()<<" "<<partt->getEnd() );
+	if(this->mStepCnt%50==49) { mpParticles->cleanup(); }
+  for(vector<ParticleObject>::iterator pit= mpParticles->getParticlesBegin();
+      pit!= mpParticles->getParticlesEnd(); pit++) {
+    //errMsg("PIT"," pit"<<(*pit)->getId()<<" pos:"<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel()  );
+    //errMsg("PIT"," pit pos:"<< (*pit).getPos()<<" vel:"<< (*pit).getVel()<<" status:"<<convertFlags2String((*pit).getFlags()) <<" " <<mpParticles->getStart()<<" "<<mpParticles->getEnd() );
 		//int flag = (*pit).getFlags();
     if( (*pit).getActive()==false ) continue;
     int i,j,k;
@@ -418,25 +402,23 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 		i= (int)(pos[0]+0.5);
 		j= (int)(pos[1]+0.5);
 		k= (int)(pos[2]+0.5);
-		if(LBMDIM==2) {
-			k = 0;
-		}
+		if(LBMDIM==2) { k = 0; }
 
 		// only testdata handling, all for sws
 #if LBM_INCLUDE_TESTSOLVERS==1
-		if(mUseTestdata) { 
-			if(mpTest->mDebugvalue1>0.0){ 
-				p->setStatus(PART_OUT);
-				mpTest->handleParticle(p, i,j,k); continue;
-		} }
+		if(useff) {
+			p->setStatus(PART_OUT);
+			mpTest->handleParticle(p, i,j,k); continue;
+		} 
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 
+		// FIXME , add k tests again, remove per type ones...
 		if(p->getStatus()&PART_IN) { // IN
 			if( (i<cutval)||(i>this->mSizex-1-cutval)||
 					(j<cutval)||(j>this->mSizey-1-cutval)
 					//||(k<cutval)||(k>this->mSizez-1-cutval) 
 					) {
-				if(!mUseTestdata) { DEL_PART;
+				if(!useff) { DEL_PART;
 				} else { 
 					p->setStatus(PART_OUT); 
 					/* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART;
@@ -454,7 +436,8 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 		}
 		//p->setStatus(PART_OUT);// DEBUG always out!
 
-		if(p->getType()==PART_BUBBLE) {
+		if( (p->getType()==PART_BUBBLE) ||
+		    (p->getType()==PART_TRACER) ) {
 
 			// no interpol
 			rho = vx = vy = vz = 0.0;
@@ -481,14 +464,14 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 					} else { // OUT
 						// out of bounds, deactivate...
 						// FIXME make fsgr treatment
-						p->setType( PART_FLOAT ); continue;
+						if(p->getType()==PART_BUBBLE) { P_CHANGETYPE(p, PART_FLOAT ); continue; }
 					}
 				} else {
 					// below 3d region, just rise
 				}
 			} else { // OUT
 #if LBM_INCLUDE_TESTSOLVERS==1
-				if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); }
+				if(useff) { mpTest->handleParticle(p, i,j,k); }
 				else DEL_PART;
 #else // LBM_INCLUDE_TESTSOLVERS==1
 				DEL_PART;
@@ -497,7 +480,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 			}
 
 			ntlVec3Gfx v = p->getVel(); // dampen...
-			if(mUseTestdata) {
+			if( (useff)&& (p->getType()==PART_BUBBLE) ) {
 				// test rise
 				//O vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
 
@@ -533,11 +516,35 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 				LbmFloat w = 0.99;
 				vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[mMaxRefine].gravity[2]);
 				v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2);
+			} else if(p->getType()==PART_TRACER) {
+				v = ntlVec3Gfx(vx,vy,vz);
+				if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid) ) {
+					// ok
+				} else {
+					const int lev = mMaxRefine;
+					LbmFloat nx,ny,nz, nv1,nv2;
+			//mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
+					if(RFLAG_NB(lev,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dE,dFfrac); } else nv1 = 0.0;
+					if(RFLAG_NB(lev,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dW,dFfrac); } else nv2 = 0.0;
+					nx = 0.5* (nv2-nv1);
+					if(RFLAG_NB(lev,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dN,dFfrac); } else nv1 = 0.0;
+					if(RFLAG_NB(lev,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dS,dFfrac); } else nv2 = 0.0;
+					ny = 0.5* (nv2-nv1);
+#if LBMDIM==3
+					if(RFLAG_NB(lev,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dT,dFfrac); } else nv1 = 0.0;
+					if(RFLAG_NB(lev,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dB,dFfrac); } else nv2 = 0.0;
+					nz = 0.5* (nv2-nv1);
+#else //LBMDIM==3
+					nz = 0.0;
+#endif //LBMDIM==3
+
+					v = (ntlVec3Gfx(nx,ny,nz)) * -0.025;
+				}
 			}
+
 			p->setVel( v );
-			//p->setVel( ntlVec3Gfx(vx,vy,vz) );
 			p->advanceVel();
-		  // fluid particle
+			//errMsg("PPPP"," pos"<<p->getPos()<<" "<<p->getVel() );
 		} 
 
 		// drop handling
@@ -583,13 +590,17 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 				if(k<cutval) { DEL_PART; continue; }
 				if(k<=this->mSizez-1-cutval){ 
 					//if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter)) {
-					if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty)) {
+					if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) {
 						// still ok
-					} else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
+					// shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
+					} else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid) ){
 						// FIXME make fsgr treatment
-						if(p->getLifeTime()>50) { 
-							p->setType( PART_FLOAT ); continue; 
-						} else DEL_PART;
+						//if(p->getLifeTime()>50) { 
+						P_CHANGETYPE(p, PART_FLOAT ); continue; 
+						// jitter in cell to prevent stacking when hitting a steep surface
+						LbmVec pos = p->getPos(); pos[0] += (rand()/(RAND_MAX+1.0))-0.5;
+						pos[1] += (rand()/(RAND_MAX+1.0))-0.5; p->setPos(pos);
+						//} else DEL_PART;
 					} else {
 						DEL_PART;
 						this->mNumParticlesLost++;
@@ -597,7 +608,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 				}
 			} else { // OUT
 #if LBM_INCLUDE_TESTSOLVERS==1
-				if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); }
+				if(useff) { mpTest->handleParticle(p, i,j,k); }
 				else{ DEL_PART; }
 #else // LBM_INCLUDE_TESTSOLVERS==1
 				{ DEL_PART; }
@@ -619,9 +630,10 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 				} else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
     			//errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
-					//p->setType( PART_BUBBLE ); continue;
+					//P_CHANGETYPE(p, PART_BUBBLE ); continue;
 					// currently bubbles off! DEBUG!
-					DEL_PART; // DEBUG bubbles off for siggraph
+					//DEL_PART; // DEBUG bubbles off for siggraph
+					P_CHANGETYPE(p, PART_FLOAT ); continue;
 				} else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
     			//errMsg("PIT","NEWDROP pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 					//? if(p->getLifeTime()>50) {
@@ -631,13 +643,14 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 							//(j<=cutval)||(j>=this->mSizey-1-cutval)) {
 					//} else 
 					//if(p->getLifeTime()>10) {
-					p->setType( PART_DROP ); continue;
+					P_CHANGETYPE(p, PART_DROP ); continue;
 					//} else DEL_PART;						
 					
 				}
 			} else { // OUT
 				// undecided particle outside... remove?
 				DEL_PART; 
+				//P_CHANGETYPE(p, PART_FLOAT ); continue;
 			}
 		}
 
@@ -650,11 +663,17 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 #if LBM_INCLUDE_TESTSOLVERS==1
 			// vanishing
 			prob = (rand()/(RAND_MAX+1.0));
-			if((mUseTestdata)&&(k>mpTest->mFluidHeight)) {
-				LbmFloat fac = (LbmFloat)(k-mpTest->mFluidHeight)/(LbmFloat)(10*(mLevel[mMaxRefine].lSizez-mpTest->mFluidHeight));
-				prob /= fac; //  TODO test? errMsg("T","T "<<prob<<" "<<fac);
+			// increse del prob. up to max height by given factor
+			const int fhCheckh = (int)(mpTest->mFluidHeight*1.25);
+			const LbmFloat fhProbh = 25.;
+			if((useff)&&(k>fhCheckh)) {
+				LbmFloat fac = (LbmFloat)(k-fhCheckh)/(LbmFloat)(fhProbh*(mLevel[mMaxRefine].lSizez-fhCheckh));
+				prob /= fac; 
 			}
 			if(prob<mLevel[mMaxRefine].timestep*0.1) DEL_PART;
+			// forced vanishing
+			//? if(k>this->mSizez*3/4) {	if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
+
 #else // LBM_INCLUDE_TESTSOLVERS==1
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 
@@ -663,13 +682,15 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 				//ntlVec3Gfx v = getVelocityAt(i,j,k);
 				rho = vx = vy = vz = 0.0;
 
-				//const int DEPTH_AVG=11; // TODO how much!?
-				const int DEPTH_AVG=7; // TODO how much!?
+				// define from particletracer.h
+#if MOVE_FLOATS==1
+				const int DEPTH_AVG=3; // only average interface vels
 				int ccnt=0;
-				for(int kk=1;kk<DEPTH_AVG;kk+=2) {
+				for(int kk=0;kk<DEPTH_AVG;kk+=1) {
 				//for(int kk=1;kk<DEPTH_AVG;kk+=1) {
 					if((k-kk)<1) continue;
-					if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue;
+					//if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue;
+					if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue;
 					ccnt++;
 					FORDF0{
 						LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l);
@@ -681,32 +702,44 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 					}
 				}
 				if(ccnt) {
-				vx /=(LbmFloat)(ccnt * 1.0); // half xy speed! value2
-				vy /=(LbmFloat)(ccnt * 1.0);
+				// use halved surface velocity (todo, use omega instead)
+				vx /=(LbmFloat)(ccnt * 2.0); // half xy speed! value2
+				vy /=(LbmFloat)(ccnt * 2.0);
 				vz /=(LbmFloat)(ccnt); }
-				// forced vanishing
-				if(k>this->mSizez*3/4) {	if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
+#else // MOVE_FLOATS==1
+				vx=vy=0.; //p->setVel(ntlVec3Gfx(0.) ); // static_float
+#endif // MOVE_FLOATS==1
+				vx += (rand()/(RAND_MAX+1.0))*FLOAT_JITTER-(FLOAT_JITTER*0.5);
+				vy += (rand()/(RAND_MAX+1.0))*FLOAT_JITTER-(FLOAT_JITTER*0.5);
 
+				bool delfloat = false;
 				if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
+					// in fluid cell
 					if((1) && (k<this->mSizez-3) && 
 							(
 							  TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFInter ) ||
-							  TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) )
-							 ) {
+							  TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) 
+							) ) {
 						vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
 						if(vz<0.0) vz=0.0;
-					} else DEL_PART;
+					} else delfloat=true;
+					// keep below obstacles
+					if((delfloat) && (TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFBnd )) ) {
+						delfloat=false; vz=0.0;
+					}
 				} else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
 					// keep in interface , one grid cell offset is added in part. gen
-				} 
-				// check if above inter, remove otherwise
-				else if((1) && (k>2) && (
-							TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) ||
-							TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) )
-							) {
-					vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
-					if(vz>0.0) vz=0.0;
-				} else DEL_PART; // */
+				} else { //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { // shipt?, was all else before
+					// check if above inter, remove otherwise
+					if((1) && (k>2) && (
+								TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) ||
+								TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) 
+								) ) {
+						vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
+						if(vz>0.0) vz=0.0;
+					} else delfloat=true; // */
+				}
+				if(delfloat) DEL_PART;
 				/*
 				// move down from empty
 				else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
@@ -715,20 +748,38 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 					//DEL_PART; // ????
 				} else  {	 DEL_PART; } // */
 				//vz = 0.0; // DEBUG
-				ntlVec3Gfx v(vx,vy,vz);
-				p->setVel( vec2G(v) ); //?
+
+				p->setVel( vec2G( ntlVec3Gfx(vx,vy,vz) ) ); //?
 				//p->setVel( vec2G(v)*0.75 + p->getVel()*0.25 ); //?
 				p->advanceVel();
 				//errMsg("PIT","IN pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 			} else {
 #if LBM_INCLUDE_TESTSOLVERS==1
-				if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); }
+				if(useff) { mpTest->handleParticle(p, i,j,k); }
 				else DEL_PART; 
 #else // LBM_INCLUDE_TESTSOLVERS==1
 				DEL_PART; 
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 				//errMsg("PIT","OUT pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 			}
+				
+			// additional bnd jitter
+			if((1) && (useff) && (p->getLifeTime()<3)) {
+				// use half butoff border 1/8
+				int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5);
+				if(maxdw<3) maxdw=3;
+#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
+#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) 
+				//if(ABS(i-(               cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
+				//if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
+				if((j>=0)&&(j<=this->mSizey-1)) {
+					if(ABS(i-(               cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(               cutval))), 0.,0.); }
+					if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); }
+				}
+//#undef FLOAT_JITTER_BND
+#undef FLOAT_JITTBNDRAND
+				//if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval)
+			}
 		} 
 		
 		// unknown particle type	
@@ -737,12 +788,13 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) {
 		}
   }
 	myTime_t parttend = getTime(); 
-	debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<" "<<partt->getNumParticles() , 10 );
+	debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<" "<<mpParticles->getNumParticles() , 10 );
 }
 
-void LbmFsgrSolver::notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename) {
+void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
 	// debug - raw dump of ffrac values
 	if(getenv("ELBEEM_RAWDEBUGDUMP")) {
+		int workSet = mLevel[mMaxRefine].setCurr;
 		std::ostringstream name;
 		//name <<"fill_" << this->mStepCnt <<".dump";
 		name << outfilename<< frameNrStr <<".dump";
@@ -752,7 +804,9 @@ void LbmFsgrSolver::notifySolverOfDump(int frameNr,char *frameNrStr,string outfi
 			for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
 				for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
 					for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
-						float val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+						float val = 0.;
+						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
 						//fwrite( &val, sizeof(val), 1, file); // binary
 						fprintf(file, "%f ",val); // text
 						//errMsg("W", PRINT_IJK<<" val:"<<val);