This version now includes the fluid control sources, however the Blender

interface for it is still missing. Right now there is only a simple hard coded
example, that moves a single control particle with strong attraction
and velocity forces through the domain.

I added more detailed information about the control code to the wiki
http://wiki.blender.org/index.php/SoCFluidDevelDoc#The_Fluid-Control_Branch ,
together with some thoughts on how a Blender integration could be done.

17 files changed, 3152 insertions, 59 deletions

diff --git a/intern/elbeem/CMakeLists.txt b/intern/elbeem/CMakeLists.txt
index 48d25901499..c0484846c6b 100644
--- a/intern/elbeem/CMakeLists.txt
+++ b/intern/elbeem/CMakeLists.txt
@@ -31,7 +31,7 @@ SET(INC ${PNG_INC} ${ZLIB_INC} ${SDL_INC})
 
 FILE(GLOB SRC intern/*.cpp)
 
-ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1)
+ADD_DEFINITIONS(-DNOGUI -DELBEEM_BLENDER=1 -DLBM_INCLUDE_CONTROL=1)
 IF(WINDOWS)
     ADD_DEFINITIONS(-DUSE_MSVC6FIXES)
 ENDIF(WINDOWS)
diff --git a/intern/elbeem/SConscript b/intern/elbeem/SConscript
index bdcb0507987..cea718ee737 100644
--- a/intern/elbeem/SConscript
+++ b/intern/elbeem/SConscript
@@ -5,7 +5,7 @@ Import('env')
 
 sources = env.Glob('intern/*.cpp')
 
-defs = ' NOGUI ELBEEM_BLENDER=1'
+defs = 'NOGUI ELBEEM_BLENDER=1 LBM_INCLUDE_CONTROL=1'
 
 if env['WITH_BF_OPENMP'] == 1:
     defs += ' PARALLEL'
diff --git a/intern/elbeem/intern/controlparticles.cpp b/intern/elbeem/intern/controlparticles.cpp
new file mode 100644
index 00000000000..cab64d7a257
--- /dev/null
+++ b/intern/elbeem/intern/controlparticles.cpp
@@ -0,0 +1,1320 @@
+// --------------------------------------------------------------------------
+//
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the 
+// GNU General Public License. See the file COPYING for details.  
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+// implementation of control particle handling
+//
+// --------------------------------------------------------------------------
+
+// indicator for LBM inclusion
+#include "ntl_geometrymodel.h"
+#include "ntl_world.h"
+#include "solver_class.h"
+#include "controlparticles.h"
+#include "mvmcoords.h"
+#include <zlib.h>
+
+#ifndef sqrtf
+#define sqrtf sqrt
+#endif
+
+// brute force circle test init in initTimeArray
+// replaced by mDebugInit
+//#define CP_FORCECIRCLEINIT 0
+
+
+void ControlParticles::initBlenderTest() {
+	mPartSets.clear();
+
+	ControlParticleSet cps;
+	mPartSets.push_back(cps);
+	int setCnt = mPartSets.size()-1;
+	ControlParticle p; 
+
+	// set for time zero
+	mPartSets[setCnt].time = 0.;
+
+	// add single particle 
+	p.reset();
+	p.pos = LbmVec(0.5, 0.5, -0.5);
+	mPartSets[setCnt].particles.push_back(p);
+
+	// add second set for animation
+	mPartSets.push_back(cps);
+	setCnt = mPartSets.size()-1;
+	mPartSets[setCnt].time = 0.15;
+
+	// insert new position
+	p.reset();
+	p.pos = LbmVec(-0.5, -0.5, 0.5);
+	mPartSets[setCnt].particles.push_back(p);
+
+	// applyTrafos();
+	initTime(0. , 1.);
+}
+
+
+// init all zero / defaults for a single particle
+void ControlParticle::reset() {
+	pos = LbmVec(0.,0.,0.);
+	vel = LbmVec(0.,0.,0.);
+	influence = 1.;
+	size = 1.;
+#ifndef LBMDIM
+#ifdef MAIN_2D
+	rotaxis = LbmVec(0.,1.,0.); // SPH xz
+#else // MAIN_2D
+	// 3d - roate in xy plane, vortex
+	rotaxis = LbmVec(0.,0.,1.);
+	// 3d - rotate for wave
+	//rotaxis = LbmVec(0.,1.,0.);
+#endif // MAIN_2D
+#else // LBMDIM
+	rotaxis = LbmVec(0.,1.,0.); // LBM xy , is swapped afterwards
+#endif // LBMDIM
+
+	density = 0.;
+	densityWeight = 0.;
+	avgVelAcc = avgVel = LbmVec(0.);
+	avgVelWeight = 0.;
+}
+
+
+// default preset/empty init
+ControlParticles::ControlParticles() :
+	_influenceTangential(0.f),
+	_influenceAttraction(0.f),
+	_influenceVelocity(0.f),
+	_influenceMaxdist(0.f),
+	_radiusAtt(1.0f),
+	_radiusVel(1.0f),
+	_radiusMinMaxd(2.0f),
+	_radiusMaxd(3.0f),
+	_currTime(-1.0), _currTimestep(1.),
+	_initTimeScale(1.), 
+	_initPartOffset(0.), _initPartScale(1.),
+	_initLastPartOffset(0.), _initLastPartScale(1.),
+	_initMirror(""),
+	_fluidSpacing(1.), _kernelWeight(-1.),
+	_charLength(1.), _charLengthInv(1.),
+	mvCPSStart(-10000.), mvCPSEnd(10000.),
+	mCPSWidth(0.1), mCPSTimestep(0.05),
+	mCPSTimeStart(0.), mCPSTimeEnd(0.5), mCPSWeightFac(1.),
+	mDebugInit(0)
+{
+	_radiusAtt = 0.15f;
+	_radiusVel = 0.15f;
+	_radiusMinMaxd = 0.16f;
+	_radiusMaxd = 0.3;
+
+	_influenceAttraction = 0.f;
+	_influenceTangential = 0.f;
+	_influenceVelocity = 0.f;
+	// 3d tests */
+}
+
+
+ 
+ControlParticles::~ControlParticles() {
+	// nothing to do...
+}
+
+LbmFloat ControlParticles::getControlTimStart() {
+	if(mPartSets.size()>0) { return mPartSets[0].time; }
+	return -1000.;
+}
+LbmFloat ControlParticles::getControlTimEnd() {
+	if(mPartSets.size()>0) { return mPartSets[mPartSets.size()-1].time; }
+	return -1000.;
+}
+
+// calculate for delta t
+void ControlParticles::setInfluenceVelocity(LbmFloat set, LbmFloat dt) {
+	const LbmFloat dtInter = 0.01;
+	LbmFloat facFv = 1.-set; //cparts->getInfluenceVelocity();
+	// mLevel[mMaxRefine].timestep
+	LbmFloat facNv = (LbmFloat)( 1.-pow( (double)facFv, (double)(dt/dtInter)) );
+	//errMsg("vwcalc","ts:"<<dt<< " its:"<<(dt/dtInter) <<" fv"<<facFv<<" nv"<<facNv<<" test:"<< pow( (double)(1.-facNv),(double)(dtInter/dt))	);
+	_influenceVelocity = facNv;
+}
+
+int ControlParticles::initExampleSet()
+{
+	// unused
+}
+
+int ControlParticles::getTotalSize()
+{
+	int s=0;
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		s+= mPartSets[i].particles.size();
+	}
+	return s;
+}
+
+// --------------------------------------------------------------------------
+// load positions & timing from text file
+// WARNING - make sure file has unix format, no win/dos linefeeds...
+#define LINE_LEN 100
+int ControlParticles::initFromTextFile(string filename)
+{
+	const bool debugRead = false;
+	char line[LINE_LEN];
+	line[LINE_LEN-1] = '\0';
+	mPartSets.clear();
+	if(filename.size()<2) return 0;
+
+	// HACK , use "cparts" suffix as old
+	// e.g. "cpart2" as new
+	if(filename[ filename.size()-1 ]=='s') {
+		return initFromTextFileOld(filename);
+	}
+
+	FILE *infile = fopen(filename.c_str(), "r");
+	if(!infile) {
+		errMsg("ControlParticles::initFromTextFile","unable to open '"<<filename<<"' " );
+		// try to open as gz sequence
+		if(initFromBinaryFile(filename)) { return 1; }
+		// try mesh MVCM generation
+		if(initFromMVCMesh(filename)) { return 1; }
+		// failed...
+		return 0;
+	}
+
+	int haveNo = false;
+	int haveScale = false;
+	int haveTime = false;
+	int noParts = -1;
+	int partCnt = 0;
+	int setCnt = 0;
+	//ControlParticle p; p.reset();
+	// scale times by constant factor while reading
+	LbmFloat timeScale= 1.0;
+	int lineCnt = 0;
+	bool abortParse = false;
+#define LASTCP mPartSets[setCnt].particles[ mPartSets[setCnt].particles.size()-1 ]
+
+	while( (!feof(infile)) && (!abortParse)) {
+		lineCnt++;
+		fgets(line, LINE_LEN, infile);
+
+		//if(debugRead) printf("\nDEBUG%d r '%s'\n",lineCnt, line);
+		if(!line) continue;
+		int len = (int)strlen(line);
+
+		// skip empty lines and comments (#,//)
+		if(len<1) continue;
+		if( (line[0]=='#') || (line[0]=='\n') ) continue;
+		if((len>1) && (line[0]=='/' && line[1]=='/')) continue;
+
+		// debug remove newline
+		if((len>=1)&&(line[len-1]=='\n')) line[len-1]='\0';
+
+		switch(line[0]) {
+
+		case 'N': { // total number of particles, more for debugging...
+			noParts = atoi(line+2);
+			if(noParts<=0) {
+				errMsg("ControlParticles::initFromTextFile","file '"<<filename<<"' - invalid no of particles "<<noParts);
+				mPartSets.clear(); fclose(infile); return 0;
+			}
+			if(debugRead) printf("CPDEBUG%d no parts '%d'\n",lineCnt, noParts );
+			haveNo = true;
+			} break;
+
+		case 'T': { // global time scale
+			timeScale *= (LbmFloat)atof(line+2);
+			if(debugRead) printf("ControlParticles::initFromTextFile - line %d , set timescale '%f', org %f\n",lineCnt, timeScale , _initTimeScale);
+			if(timeScale==0.) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error: timescale = 0.! reseting to 1 ...\n",lineCnt); timeScale=1.; }
+			haveScale = true;
+			} break;
+
+		case 'I': { // influence settings, overrides others as of now...
+			float val = (LbmFloat)atof(line+3);
+			const char *setvar = "[invalid]";
+			switch(line[1]) {
+				//case 'f': { _influenceFalloff = val; setvar = "falloff"; } break;
+				case 't': { _influenceTangential = val; setvar = "tangential"; } break;
+				case 'a': { _influenceAttraction = val; setvar = "attraction"; } break;
+				case 'v': { _influenceVelocity = val; setvar = "velocity"; } break;
+				case 'm': { _influenceMaxdist = val; setvar = "maxdist"; } break;
+				default: 
+					fprintf(stdout,"ControlParticles::initFromTextFile (%s) - line %d , invalid influence setting %c, %f\n",filename.c_str() ,lineCnt, line[1], val);
+			}
+			if(debugRead) printf("CPDEBUG%d set influence '%s'=%f \n",lineCnt, setvar, val);
+			} break;
+
+		case 'R': { // radius settings, overrides others as of now...
+			float val = (LbmFloat)atof(line+3);
+			const char *setvar = "[invalid]";
+			switch(line[1]) {
+				case 'a': { _radiusAtt = val; setvar = "r_attraction"; } break;
+				case 'v': { _radiusVel = val; setvar = "r_velocity"; } break;
+				case 'm': { _radiusMaxd = val; setvar = "r_maxdist"; } break;
+				default: 
+					fprintf(stdout,"ControlParticles::initFromTextFile (%s) - line %d , invalid influence setting %c, %f\n",filename.c_str() ,lineCnt, line[1], val);
+			}
+			if(debugRead) printf("CPDEBUG%d set influence '%s'=%f \n",lineCnt, setvar, val);
+			} break;
+
+		case 'S': { // new particle set at time T
+			ControlParticleSet cps;
+			mPartSets.push_back(cps);
+			setCnt = (int)mPartSets.size()-1;
+
+			LbmFloat val = (LbmFloat)atof(line+2);
+			mPartSets[setCnt].time = val * timeScale;
+			if(debugRead) printf("CPDEBUG%d new set, time '%f', %d\n",lineCnt, mPartSets[setCnt].time, setCnt );
+			haveTime = true;
+			partCnt = -1;
+			} break;
+
+		case 'P':   // new particle with pos
+		case 'n': { // new particle without pos
+				if((!haveTime)||(setCnt<0)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error: set missing!\n",lineCnt); abortParse=true; break; }
+				partCnt++;
+				if(partCnt>=noParts) {
+					if(debugRead) printf("CPDEBUG%d partset done \n",lineCnt);
+					haveTime = false;
+				} else {
+					ControlParticle p; p.reset();
+					mPartSets[setCnt].particles.push_back(p);
+				}
+			} 
+			// only new part, or new with pos?
+			if(line[0] == 'n') break;
+
+		// particle properties
+
+		case 'p': { // new particle set at time T
+			if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|p: particle missing!\n",lineCnt); abortParse=true; break; }
+			float px=0.,py=0.,pz=0.;
+			if( sscanf(line+2,"%f %f %f",&px,&py,&pz) != 3) {
+				fprintf(stdout,"CPDEBUG%d, unable to parse position!\n",lineCnt); abortParse=true; break; 
+			}
+			if(!(finite(px)&&finite(py)&&finite(pz))) { px=py=pz=0.; }
+			LASTCP.pos[0] = px;
+			LASTCP.pos[1] = py;
+			LASTCP.pos[2] = pz; 
+			if(debugRead) printf("CPDEBUG%d part%d,%d: position %f,%f,%f \n",lineCnt,setCnt,partCnt, px,py,pz);
+			} break;
+
+		case 's': { // particle size
+			if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|s: particle missing!\n",lineCnt); abortParse=true; break; }
+			float ps=1.;
+			if( sscanf(line+2,"%f",&ps) != 1) {
+				fprintf(stdout,"CPDEBUG%d, unable to parse size!\n",lineCnt); abortParse=true; break; 
+			}
+			if(!(finite(ps))) { ps=0.; }
+			LASTCP.size = ps;
+			if(debugRead) printf("CPDEBUG%d part%d,%d: size %f \n",lineCnt,setCnt,partCnt, ps);
+			} break;
+
+		case 'i': { // particle influence
+			if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|i: particle missing!\n",lineCnt); abortParse=true; break; }
+			float pinf=1.;
+			if( sscanf(line+2,"%f",&pinf) != 1) {
+				fprintf(stdout,"CPDEBUG%d, unable to parse size!\n",lineCnt); abortParse=true; break; 
+			}
+			if(!(finite(pinf))) { pinf=0.; }
+			LASTCP.influence = pinf;
+			if(debugRead) printf("CPDEBUG%d part%d,%d: influence %f \n",lineCnt,setCnt,partCnt, pinf);
+			} break;
+
+		case 'a': { // rotation axis
+			if((!haveTime)||(setCnt<0)||(mPartSets[setCnt].particles.size()<1)) { fprintf(stdout,"ControlParticles::initFromTextFile - line %d ,error|a: particle missing!\n",lineCnt); abortParse=true; break; }
+			float px=0.,py=0.,pz=0.;
+			if( sscanf(line+2,"%f %f %f",&px,&py,&pz) != 3) {
+				fprintf(stdout,"CPDEBUG%d, unable to parse rotaxis!\n",lineCnt); abortParse=true; break; 
+			}
+			if(!(finite(px)&&finite(py)&&finite(pz))) { px=py=pz=0.; }
+			LASTCP.rotaxis[0] = px;
+			LASTCP.rotaxis[1] = py;
+			LASTCP.rotaxis[2] = pz; 
+			if(debugRead) printf("CPDEBUG%d part%d,%d: rotaxis %f,%f,%f \n",lineCnt,setCnt,partCnt, px,py,pz);
+			} break;
+
+
+		default:
+			if(debugRead) printf("CPDEBUG%d ignored: '%s'\n",lineCnt, line );
+			break;
+		}
+	}
+	if(debugRead && abortParse) printf("CPDEBUG aborted parsing after set... %d\n",(int)mPartSets.size() );
+
+	// sanity check
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		if( (int)mPartSets[i].particles.size()!=noParts) {
+			fprintf(stdout,"ControlParticles::initFromTextFile (%s) - invalid no of particles in set %d, is:%d, shouldbe:%d \n",filename.c_str() ,i,(int)mPartSets[i].particles.size(), noParts);
+			mPartSets.clear();
+			fclose(infile);
+			return 0;
+		}
+	}
+
+	// print stats
+	printf("ControlParticles::initFromTextFile (%s): Read %d sets, each %d particles\n",filename.c_str() ,
+			(int)mPartSets.size(), noParts );
+	if(mPartSets.size()>0) {
+		printf("ControlParticles::initFromTextFile (%s): Time: %f,%f\n",filename.c_str() ,mPartSets[0].time, mPartSets[mPartSets.size()-1].time );
+	}
+	
+	// done...
+	fclose(infile);
+	applyTrafos();
+	return 1;
+}
+
+
+int ControlParticles::initFromTextFileOld(string filename)
+{
+	const bool debugRead = false;
+	char line[LINE_LEN];
+	line[LINE_LEN-1] = '\0';
+	mPartSets.clear();
+	if(filename.size()<1) return 0;
+
+	FILE *infile = fopen(filename.c_str(), "r");
+	if(!infile) {
+		fprintf(stdout,"ControlParticles::initFromTextFileOld - unable to open '%s'\n",filename.c_str() );
+		return 0;
+	}
+
+	int haveNo = false;
+	int haveScale = false;
+	int haveTime = false;
+	int noParts = -1;
+	int coordCnt = 0;
+	int partCnt = 0;
+	int setCnt = 0;
+	ControlParticle p; p.reset();
+	// scale times by constant factor while reading
+	LbmFloat timeScale= 1.0;
+	int lineCnt = 0;
+
+	while(!feof(infile)) {
+		lineCnt++;
+		fgets(line, LINE_LEN, infile);
+
+		if(debugRead) printf("\nDEBUG%d r '%s'\n",lineCnt, line);
+
+		if(!line) continue;
+		int len = (int)strlen(line);
+
+		// skip empty lines and comments (#,//)
+		if(len<1) continue;
+		if( (line[0]=='#') || (line[0]=='\n') ) continue;
+		if((len>1) && (line[0]=='/' && line[1]=='/')) continue;
+
+		// debug remove newline
+		if((len>=1)&&(line[len-1]=='\n')) line[len-1]='\0';
+
+		// first read no. of particles
+		if(!haveNo) {
+			noParts = atoi(line);
+			if(noParts<=0) {
+				fprintf(stdout,"ControlParticles::initFromTextFileOld - invalid no of particles %d\n",noParts);
+				mPartSets.clear();
+				fclose(infile);
+				return 0;
+			}
+			if(debugRead) printf("DEBUG%d noparts '%d'\n",lineCnt, noParts );
+			haveNo = true;
+		} 
+
+		// then read time scale
+		else if(!haveScale) {
+			timeScale *= (LbmFloat)atof(line);
+			if(debugRead) printf("DEBUG%d tsc '%f', org %f\n",lineCnt, timeScale , _initTimeScale);
+			haveScale = true;
+		} 
+
+		// then get set time
+		else if(!haveTime) {
+			ControlParticleSet cps;
+			mPartSets.push_back(cps);
+			setCnt = (int)mPartSets.size()-1;
+
+			LbmFloat val = (LbmFloat)atof(line);
+			mPartSets[setCnt].time = val * timeScale;
+			if(debugRead) printf("DEBUG%d time '%f', %d\n",lineCnt, mPartSets[setCnt].time, setCnt );
+			haveTime = true;
+		}
+
+		// default read all parts
+		else {
+			LbmFloat val = (LbmFloat)atof(line);
+			if(debugRead) printf("DEBUG: l%d s%d,particle%d '%f' %d,%d/%d\n",lineCnt,(int)mPartSets.size(),(int)mPartSets[setCnt].particles.size(), val ,coordCnt,partCnt,noParts);
+			p.pos[coordCnt] = val;
+			coordCnt++;
+			if(coordCnt>=3) {
+				mPartSets[setCnt].particles.push_back(p);
+				p.reset();
+				coordCnt=0;
+				partCnt++;
+			}
+			if(partCnt>=noParts) {
+				partCnt = 0;
+				haveTime = false;
+			}
+			//if(debugRead) printf("DEBUG%d par2 %d,%d/%d\n",lineCnt, coordCnt,partCnt,noParts);
+		}
+		//read pos, vel ...
+	}
+
+	// sanity check
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		if( (int)mPartSets[i].particles.size()!=noParts) {
+			fprintf(stdout,"ControlParticles::initFromTextFileOld - invalid no of particles in set %d, is:%d, shouldbe:%d \n",i,(int)mPartSets[i].particles.size(), noParts);
+			mPartSets.clear();
+			fclose(infile);
+			return 0;
+		}
+	}
+	// print stats
+	printf("ControlParticles::initFromTextFileOld: Read %d sets, each %d particles\n",
+			(int)mPartSets.size(), noParts );
+	if(mPartSets.size()>0) {
+		printf("ControlParticles::initFromTextFileOld: Time: %f,%f\n",mPartSets[0].time, mPartSets[mPartSets.size()-1].time );
+	}
+	
+	// done...
+	fclose(infile);
+	applyTrafos();
+	return 1;
+}
+
+// load positions & timing from gzipped binary file
+int ControlParticles::initFromBinaryFile(string filename) {
+	mPartSets.clear();
+	if(filename.size()<1) return 0;
+	int fileNotFound=0;
+	int fileFound=0;
+	char ofile[256];
+
+	for(int set=0; ((set<10000)&&(fileNotFound<10)); set++) {
+		snprintf(ofile,256,"%s%04d.gz",filename.c_str(),set);
+		//errMsg("ControlParticle::initFromBinaryFile","set"<<set<<" notf"<<fileNotFound<<" ff"<<fileFound);
+
+		gzFile gzf;
+		gzf = gzopen(ofile, "rb");
+		if (!gzf) {
+			//errMsg("ControlParticles::initFromBinaryFile","Unable to open file for reading '"<<ofile<<"' "); 
+			fileNotFound++;
+			continue;
+		}
+		fileNotFound=0;
+		fileFound++;
+
+		ControlParticleSet cps;
+		mPartSets.push_back(cps);
+		int setCnt = (int)mPartSets.size()-1;
+		//LbmFloat val = (LbmFloat)atof(line+2);
+		mPartSets[setCnt].time = (gfxReal)set;
+
+		int totpart = 0;
+		gzread(gzf, &totpart, sizeof(totpart));
+
+		for(int a=0; a<totpart; a++) {
+			int ptype=0;
+			float psize=0.0;
+			ntlVec3Gfx ppos,pvel;
+			gzread(gzf, &ptype, sizeof( ptype )); 
+			gzread(gzf, &psize, sizeof( float )); 
+
+			for(int j=0; j<3; j++) { gzread(gzf, &ppos[j], sizeof( float )); }
+			for(int j=0; j<3; j++) { gzread(gzf, &pvel[j], sizeof( float )); }
+
+			ControlParticle p; 
+			p.reset();
+			p.pos = vec2L(ppos);
+			mPartSets[setCnt].particles.push_back(p);
+		} 
+
+		gzclose(gzf);
+		//errMsg("ControlParticle::initFromBinaryFile","Read set "<<ofile<<", #"<<mPartSets[setCnt].particles.size() ); // DEBUG
+	} // sets
+
+	if(fileFound==0) return 0;
+	applyTrafos();
+	return 1;
+}
+
+int globCPIProblems =0;
+bool ControlParticles::checkPointInside(ntlTree *tree, ntlVec3Gfx org, gfxReal &distance) {
+	// warning - stripped down version of geoInitCheckPointInside
+	const int globGeoInitDebug = 0;
+	const int  flags = FGI_FLUID;
+	org += ntlVec3Gfx(0.0001);
+	ntlVec3Gfx dir = ntlVec3Gfx(1.0, 0.0, 0.0);
+	int OId = -1;
+	ntlRay ray(org, dir, 0, 1.0, NULL);
+	bool done = false;
+	bool inside = false;
+	int mGiObjInside = 0; 
+	LbmFloat mGiObjDistance = -1.0; 
+	LbmFloat giObjFirstHistSide = 0; 
+	
+	// if not inside, return distance to first hit
+	gfxReal firstHit=-1.0;
+	int     firstOId = -1;
+	if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org);
+
+	while(!done) {
+		// find first inside intersection
+		ntlTriangle *triIns = NULL;
+		distance = -1.0;
+		ntlVec3Gfx normal(0.0);
+		tree->intersectX(ray,distance,normal, triIns, flags, true);
+		if(triIns) {
+			ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
+			LbmFloat orientation = dot(normal, dir);
+			OId = triIns->getObjectId();
+			if(orientation<=0.0) {
+				// outside hit
+				normal *= -1.0;
+				mGiObjInside++;
+				if(giObjFirstHistSide==0) giObjFirstHistSide = 1;
+				if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
+			} else {
+				// inside hit
+				mGiObjInside++;
+				if(mGiObjDistance<0.0) mGiObjDistance = distance;
+				if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
+				if(giObjFirstHistSide==0) giObjFirstHistSide = -1;
+			}
+			norg += normal * getVecEpsilon();
+			ray = ntlRay(norg, dir, 0, 1.0, NULL);
+			// remember first hit distance, in case we're not 
+			// inside anything
+			if(firstHit<0.0) {
+				firstHit = distance;
+				firstOId = OId;
+			}
+		} else {
+			// no more intersections... return false
+			done = true;
+		}
+	}
+
+	distance = -1.0;
+	if(mGiObjInside>0) {
+		bool mess = false;
+		if((mGiObjInside%2)==1) {
+			if(giObjFirstHistSide != -1) mess=true;
+		} else {
+			if(giObjFirstHistSide !=  1) mess=true;
+		}
+		if(mess) {
+			// ?
+			//errMsg("IIIproblem","At "<<org<<" obj  inside:"<<mGiObjInside<<" firstside:"<<giObjFirstHistSide );
+			globCPIProblems++;
+			mGiObjInside++; // believe first hit side...
+		}
+	}
+
+	if(globGeoInitDebug) errMsg("CHIII"," ins="<<mGiObjInside<<" t"<<mGiObjDistance<<" d"<<distance);
+	if(((mGiObjInside%2)==1)&&(mGiObjDistance>0.0)) {
+		if(  (distance<0.0)                             || // first intersection -> good
+				((distance>0.0)&&(distance>mGiObjDistance)) // more than one intersection -> use closest one
+			) {						
+			distance = mGiObjDistance;
+			OId = 0;
+			inside = true;
+		} 
+	}
+
+	if(!inside) {
+		distance = firstHit;
+		OId = firstOId;
+	}
+	if(globGeoInitDebug) errMsg("CHIII","ins"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
+
+	return inside;
+}
+int ControlParticles::initFromMVCMesh(string filename) {
+	myTime_t mvmstart = getTime(); 
+	ntlGeometryObjModel *model = new ntlGeometryObjModel();
+	int gid=1;
+	char infile[256];
+	vector<ntlTriangle> triangles;
+	vector<ntlVec3Gfx> vertices;
+	vector<ntlVec3Gfx> normals;
+	snprintf(infile,256,"%s.bobj.gz", filename.c_str() );
+	model->loadBobjModel(string(infile));
+	model->setLoaded(true);
+	model->setGeoInitId(gid);
+	model->setGeoInitType(FGI_FLUID);
+	debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"infile:"<<string(infile) ,4);
+
+	//getTriangles(double t,  vector<ntlTriangle> *triangles, vector<ntlVec3Gfx> *vertices, vector<ntlVec3Gfx> *normals, int objectId );
+	model->getTriangles(mCPSTimeStart, &triangles, &vertices, &normals, 1 ); 
+	debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG," tris:"<<triangles.size()<<" verts:"<<vertices.size()<<" norms:"<<normals.size() , 2);
+	
+	// valid mesh?
+	if(triangles.size() <= 0) {
+		return 0;
+	}
+
+	ntlRenderGlobals *glob = new ntlRenderGlobals;
+ 	ntlScene *genscene = new ntlScene( glob, false );
+	genscene->addGeoClass(model);
+	genscene->addGeoObject(model);
+	genscene->buildScene(0., false);
+	char treeFlag = (1<<(4+gid));
+
+	ntlTree *tree = new ntlTree( 
+			15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
+			genscene, treeFlag );
+
+	// TODO? use params
+	ntlVec3Gfx start,end;
+	model->getExtends(start,end);
+
+	LbmFloat width = mCPSWidth;
+	if(width<=LBM_EPSILON) { errMsg("ControlParticles::initFromMVMCMesh","Invalid mCPSWidth! "<<mCPSWidth); width=mCPSWidth=0.1; }
+	ntlVec3Gfx org = start+ntlVec3Gfx(width*0.5);
+	gfxReal distance = -1.;
+	vector<ntlVec3Gfx> inspos;
+	int approxmax = (int)( ((end[0]-start[0])/width)*((end[1]-start[1])/width)*((end[2]-start[2])/width) );
+
+	debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"start"<<start<<" end"<<end<<" w="<<width<<" maxp:"<<approxmax, 5);
+	while(org[2]<end[2]) {
+		while(org[1]<end[1]) {
+			while(org[0]<end[0]) {
+				if(checkPointInside(tree, org, distance)) {
+					inspos.push_back(org);
+					//inspos.push_back(org+ntlVec3Gfx(width));
+					//inspos.push_back(start+end*0.5);
+				}
+				// TODO optimize, use distance
+				org[0] += width;
+			}
+			org[1] += width;
+			org[0] = start[0];
+		}
+		org[2] += width;
+		org[1] = start[1];
+	}
+	debMsgStd("ControlParticles::initFromMVMCMesh",DM_MSG,"points: "<<inspos.size()<<" initproblems: "<<globCPIProblems,5 );
+
+	MeanValueMeshCoords mvm;
+	mvm.calculateMVMCs(vertices,triangles, inspos, mCPSWeightFac);
+	vector<ntlVec3Gfx> ninspos;
+	mvm.transfer(vertices, ninspos);
+
+	// init first set, check dist
+	ControlParticleSet firstcps; //T
+	mPartSets.push_back(firstcps);
+	mPartSets[mPartSets.size()-1].time = (gfxReal)0.;
+	vector<bool> useCP;
+	bool debugPos=false;
+
+	for(int i=0; i<(int)inspos.size(); i++) {
+		ControlParticle p; p.reset();
+		p.pos = vec2L(inspos[i]);
+		//errMsg("COMP "," "<<inspos[i]<<" vs "<<ninspos[i] );
+		double cpdist = norm(inspos[i]-ninspos[i]);
+		bool usecpv = true;
+		if(debugPos) errMsg("COMP "," "<<cpdist<<usecpv);
+
+		mPartSets[mPartSets.size()-1].particles.push_back(p);
+		useCP.push_back(usecpv);
+	}
+
+	// init further sets, temporal mesh sampling
+	double tsampling = mCPSTimestep;
+	int totcnt = (int)( (mCPSTimeEnd-mCPSTimeStart)/tsampling ), tcnt=0;
+	for(double t=mCPSTimeStart+tsampling; ((t<mCPSTimeEnd) && (ninspos.size()>0.)); t+=tsampling) {
+		ControlParticleSet nextcps; //T
+		mPartSets.push_back(nextcps);
+		mPartSets[mPartSets.size()-1].time = (gfxReal)t;
+
+		vertices.clear(); triangles.clear(); normals.clear();
+		model->getTriangles(t, &triangles, &vertices, &normals, 1 );
+		mvm.transfer(vertices, ninspos);
+		if(tcnt%(totcnt/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Transferring animation, frame: "<<tcnt<<"/"<<totcnt,5 );
+		tcnt++;
+		for(int i=0; i<(int)ninspos.size(); i++) {
+			if(debugPos) errMsg("COMP "," "<<norm(inspos[i]-ninspos[i]) );
+			if(useCP[i]) {
+				ControlParticle p; p.reset();
+				p.pos = vec2L(ninspos[i]);
+				mPartSets[mPartSets.size()-1].particles.push_back(p);
+			}
+		}
+	}
+
+	applyTrafos();
+
+	myTime_t mvmend = getTime(); 
+	debMsgStd("ControlParticle::initFromMVMCMesh",DM_MSG,"t:"<<getTimeString(mvmend-mvmstart)<<" ",7 );
+	delete tree;
+	delete genscene;
+	delete glob;
+//exit(1); // DEBUG
+	return 1;
+}
+
+#define TRISWAP(v,a,b) { LbmFloat tmp = (v)[b]; (v)[b]=(v)[a]; (v)[a]=tmp; }
+#define TRISWAPALL(v,a,b) {  \
+			TRISWAP( (v).pos     ,a,b ); \
+			TRISWAP( (v).vel     ,a,b ); \
+			TRISWAP( (v).rotaxis ,a,b ); }
+
+// helper function for LBM 2D -> swap Y and Z components everywhere
+void ControlParticles::swapCoords(int a, int b) {
+	//return;
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		for(int j=0; j<(int)mPartSets[i].particles.size(); j++) {
+			TRISWAPALL( mPartSets[i].particles[j],a,b );
+		}
+	}
+}
+
+// helper function for LBM 2D -> mirror time
+void ControlParticles::mirrorTime() {
+	LbmFloat maxtime = mPartSets[mPartSets.size()-1].time;
+	const bool debugTimeswap = false;
+	
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		mPartSets[i].time = maxtime - mPartSets[i].time;
+	}
+
+	for(int i=0; i<(int)mPartSets.size()/2; i++) {
+		ControlParticleSet cps = mPartSets[i];
+		if(debugTimeswap) errMsg("TIMESWAP", " s"<<i<<","<<mPartSets[i].time<<"  and s"<<(mPartSets.size()-1-i)<<","<< mPartSets[mPartSets.size()-1-i].time <<"  mt:"<<maxtime );
+		mPartSets[i] = mPartSets[mPartSets.size()-1-i];
+		mPartSets[mPartSets.size()-1-i] = cps;
+	}
+
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		if(debugTimeswap) errMsg("TIMESWAP", "done: s"<<i<<","<<mPartSets[i].time<<"  "<<mPartSets[i].particles.size() );
+	}
+}
+
+// apply init transformations
+void ControlParticles::applyTrafos() {
+	// apply trafos
+	for(int i=0; i<(int)mPartSets.size(); i++) {
+		mPartSets[i].time *= _initTimeScale;
+		/*for(int j=0; j<(int)mPartSets[i].particles.size(); j++) {
+			for(int k=0; k<3; k++) {
+				mPartSets[i].particles[j].pos[k] *= _initPartScale[k];
+				mPartSets[i].particles[j].pos[k] += _initPartOffset[k];
+			}
+		} now done in initarray */
+	}
+
+	// mirror coords...
+	for(int l=0; l<(int)_initMirror.length(); l++) {
+		switch(_initMirror[l]) {
+		case 'X':
+		case 'x':
+			//printf("ControlParticles::applyTrafos - mirror x\n");
+			swapCoords(1,2);
+			break;
+		case 'Y':
+		case 'y':
+			//printf("ControlParticles::applyTrafos - mirror y\n");
+			swapCoords(0,2);
+			break;
+		case 'Z':
+		case 'z':
+			//printf("ControlParticles::applyTrafos - mirror z\n");
+			swapCoords(0,1);
+			break;
+		case 'T':
+		case 't':
+			//printf("ControlParticles::applyTrafos - mirror time\n");
+			mirrorTime();
+			break;
+		case ' ':
+		case '-':
+		case '\n':
+			break;
+		default:
+			//printf("ControlParticles::applyTrafos - mirror unknown %c !?\n", _initMirror[l] );
+			break;
+		}
+	}
+
+	// reset 2d positions
+#if (CP_PROJECT2D==1) && ( defined(MAIN_2D) || LBMDIM==2 )
+	for(size_t j=0; j<mPartSets.size(); j++) 
+		for(size_t i=0; i<mPartSets[j].particles.size(); i++) {
+			// DEBUG 
+			mPartSets[j].particles[i].pos[1] = 0.f;
+		}
+#endif
+
+#if defined(LBMDIM) 
+	//? if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){ 
+		// gui control test, don swap...
+	//? } else {
+		//? swapCoords(1,2); // LBM 2D -> swap Y and Z components everywhere
+	//? }
+#endif
+
+	initTime(0.f, 0.f);
+}
+
+#undef TRISWAP
+
+// --------------------------------------------------------------------------
+// init for a given time
+void ControlParticles::initTime(LbmFloat t, LbmFloat dt) 
+{
+	//fprintf(stdout, "CPINITTIME init %f\n",t);
+	_currTime = t;
+	if(mPartSets.size()<1) return;
+
+	// init zero velocities
+	initTimeArray(t, _particles);
+
+	// calculate velocities from prev. timestep?
+	if(dt>0.) {
+		_currTimestep = dt;
+		std::vector<ControlParticle> prevparts;
+		initTimeArray(t-dt, prevparts);
+		LbmFloat invdt = 1.0/dt;
+		for(size_t j=0; j<_particles.size(); j++) {
+			ControlParticle &p = _particles[j];
+			ControlParticle &prevp = prevparts[j];
+			for(int k=0; k<3; k++) {
+				p.pos[k] *= _initPartScale[k];
+				p.pos[k] += _initPartOffset[k];
+				prevp.pos[k] *= _initLastPartScale[k];
+				prevp.pos[k] += _initLastPartOffset[k];
+			}
+			p.vel = (p.pos - prevp.pos)*invdt;
+		}
+
+		if(0) {
+			LbmVec avgvel(0.);
+			for(size_t j=0; j<_particles.size(); j++) {
+				avgvel += _particles[j].vel;
+			}
+			avgvel /= (LbmFloat)_particles.size();
+			//fprintf(stdout," AVGVEL %f,%f,%f \n",avgvel[0],avgvel[1],avgvel[2]); // DEBUG
+		}
+	}
+}
+
+// helper, init given array
+void ControlParticles::initTimeArray(LbmFloat t, std::vector<ControlParticle> &parts) {
+	if(mPartSets.size()<1) return;
+
+	if(parts.size()!=mPartSets[0].particles.size()) {
+		//fprintf(stdout,"PRES \n");
+		parts.resize(mPartSets[0].particles.size());
+		// TODO reset all?
+		for(size_t j=0; j<parts.size(); j++) {
+			parts[j].reset();
+		}
+	}
+	if(parts.size()<1) return;
+
+	// debug inits
+	if(mDebugInit==1) {
+		// hard coded circle init
+		for(size_t j=0; j<mPartSets[0].particles.size(); j++) {
+			ControlParticle p = mPartSets[0].particles[j];
+			// remember old
+			p.density = parts[j].density;
+			p.densityWeight = parts[j].densityWeight;
+			p.avgVel = parts[j].avgVel;
+			p.avgVelAcc = parts[j].avgVelAcc;
+			p.avgVelWeight = parts[j].avgVelWeight;
+			LbmVec ppos(0.); { // DEBUG
+			const float tscale=10.;
+			const float tprevo = 0.33;
+			const LbmVec toff(50,50,0);
+			const LbmVec oscale(30,30,0);
+			ppos[0] =  cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0] + toff[0];
+			ppos[1] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1] + toff[1];
+			ppos[2] =                               toff[2]; } // DEBUG
+			p.pos = ppos;
+			parts[j] = p;
+			//errMsg("ControlParticle::initTimeArray","j:"<<j<<" p:"<<parts[j].pos );
+		}
+		return;
+	}
+	else if(mDebugInit==2) {
+		// hard coded spiral init
+		const float tscale=-10.;
+		const float tprevo = 0.33;
+		LbmVec   toff(50,0,-50);
+		const LbmVec oscale(20,20,0);
+		toff[2] += 30. * t +30.;
+		for(size_t j=0; j<mPartSets[0].particles.size(); j++) {
+			ControlParticle p = mPartSets[0].particles[j];
+			// remember old
+			p.density = parts[j].density;
+			p.densityWeight = parts[j].densityWeight;
+			p.avgVel = parts[j].avgVel;
+			p.avgVelAcc = parts[j].avgVelAcc;
+			p.avgVelWeight = parts[j].avgVelWeight;
+			LbmVec ppos(0.); 
+			ppos[1] =                               toff[2]; 
+			LbmFloat zscal = (ppos[1]+100.)/200.;
+			ppos[0] =  cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0]*zscal + toff[0];
+			ppos[2] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1]*zscal + toff[1];
+			p.pos = ppos;
+			parts[j] = p;
+
+			toff[2] += 0.25;
+		}
+		return;
+	}
+
+	// use first set
+	if((t<=mPartSets[0].time)||(mPartSets.size()==1)) {
+		//fprintf(stdout,"PINI %f \n", t);
+		//parts = mPartSets[0].particles;
+		const int i=0;
+		for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+			ControlParticle p = mPartSets[i].particles[j];
+			// remember old
+			p.density = parts[j].density;
+			p.densityWeight = parts[j].densityWeight;
+			p.avgVel = parts[j].avgVel;
+			p.avgVelAcc = parts[j].avgVelAcc;
+			p.avgVelWeight = parts[j].avgVelWeight;
+			parts[j] = p;
+		}
+		return;
+	}
+
+	for(int i=0; i<(int)mPartSets.size()-1; i++) {
+		if((mPartSets[i].time<=t) && (mPartSets[i+1].time>t)) {
+			LbmFloat d = mPartSets[i+1].time-mPartSets[i].time;
+			LbmFloat f = (t-mPartSets[i].time)/d;
+			LbmFloat omf = 1.0f - f;
+	
+			for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+				ControlParticle *src1=&mPartSets[i  ].particles[j];
+				ControlParticle *src2=&mPartSets[i+1].particles[j];
+				ControlParticle &p = parts[j];
+				// do linear interpolation
+				p.pos     = src1->pos * omf     + src2->pos *f;
+				p.vel     = LbmVec(0.); // reset, calculated later on src1->vel * omf     + src2->vel *f;
+				p.rotaxis = src1->rotaxis * omf + src2->rotaxis *f;
+				p.influence = src1->influence * omf + src2->influence *f;
+				p.size    = src1->size * omf    + src2->size *f;
+				// dont modify: density, densityWeight
+			}
+		}
+	}
+
+	// after last?
+	if(t>=mPartSets[ mPartSets.size() -1 ].time) {
+		//parts = mPartSets[ mPartSets.size() -1 ].particles;
+		const int i= (int)mPartSets.size() -1;
+		for(size_t j=0; j<mPartSets[i].particles.size(); j++) {
+			ControlParticle p = mPartSets[i].particles[j];
+			// restore
+			p.density = parts[j].density;
+			p.densityWeight = parts[j].densityWeight;
+			p.avgVel = parts[j].avgVel;
+			p.avgVelAcc = parts[j].avgVelAcc;
+			p.avgVelWeight = parts[j].avgVelWeight;
+			parts[j] = p;
+		}
+	}
+}
+
+
+
+
+// --------------------------------------------------------------------------
+
+#define DEBUG_MODVEL 0
+
+// recalculate 
+void ControlParticles::calculateKernelWeight() {
+	const bool debugKernel = true;
+
+	// calculate kernel area with respect to particlesize/cellsize
+	LbmFloat kernelw = -1.;
+	LbmFloat kernelnorm = -1.;
+	LbmFloat krad = (_radiusAtt*0.75); // FIXME  use real cone approximation...?
+	//krad = (_influenceFalloff*1.);
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+	kernelw = CP_PI*krad*krad;
+	kernelnorm = 1.0 / (_fluidSpacing * _fluidSpacing);
+#else // 2D
+	kernelw = CP_PI*krad*krad*krad* (4./3.);
+	kernelnorm = 1.0 / (_fluidSpacing * _fluidSpacing * _fluidSpacing);
+#endif // MAIN_2D
+
+	if(debugKernel) debMsgStd("ControlParticles::calculateKernelWeight",DM_MSG,"kw"<<kernelw<<", norm"<<
+			kernelnorm<<", w*n="<<(kernelw*kernelnorm)<<", rad"<<krad<<", sp"<<_fluidSpacing<<"  ", 7);
+	LbmFloat kernelws = kernelw*kernelnorm;
+	_kernelWeight = kernelws;
+	if(debugKernel) debMsgStd("ControlParticles::calculateKernelWeight",DM_MSG,"influence f="<<_radiusAtt<<" t="<<
+			_influenceTangential<<" a="<<_influenceAttraction<<" v="<<_influenceVelocity<<" kweight="<<_kernelWeight, 7);
+	if(_kernelWeight<=0.) {
+		errMsg("ControlParticles::calculateKernelWeight", "invalid kernel! "<<_kernelWeight<<", resetting");
+		_kernelWeight = 1.;
+	}
+}
+
+void 
+ControlParticles::prepareControl(LbmFloat simtime, LbmFloat dt, ControlParticles *motion) {
+	debMsgStd("ControlParticle::prepareControl",DM_MSG," simtime="<<simtime<<" dt="<<dt<<" ", 5);
+
+	//fprintf(stdout,"PREPARE \n");
+	LbmFloat avgdw = 0.;
+	for(size_t i=0; i<_particles.size(); i++) {
+		ControlParticle *cp = &_particles[i];
+
+		if(this->getInfluenceAttraction()<0.) {
+			cp->density= 
+			cp->densityWeight = 1.0;
+			continue;
+		} 
+
+		// normalize by kernel 
+		//cp->densityWeight = (1.0 - (cp->density / _kernelWeight)); // store last
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+		cp->densityWeight = (1.0 - (cp->density / (_kernelWeight*cp->size*cp->size) )); // store last
+#else // 2D
+		cp->densityWeight = (1.0 - (cp->density / (_kernelWeight*cp->size*cp->size*cp->size) )); // store last
+#endif // MAIN_2D
+
+		if(i<10) debMsgStd("ControlParticle::prepareControl",DM_MSG,"kernelDebug i="<<i<<" densWei="<<cp->densityWeight<<" 1/kw"<<(1.0/_kernelWeight)<<" cpdensity="<<cp->density, 9 );
+		if(cp->densityWeight<0.) cp->densityWeight=0.;
+		if(cp->densityWeight>1.) cp->densityWeight=1.;
+		
+		avgdw += cp->densityWeight;
+		// reset for next step
+		cp->density = 0.; 
+
+		if(cp->avgVelWeight>0.) {
+		 cp->avgVel	= cp->avgVelAcc/cp->avgVelWeight; 
+		 cp->avgVelWeight	= 0.; 
+		 cp->avgVelAcc	= LbmVec(0.,0.,0.); 
+		}
+	}
+	//if(debugKernel) for(size_t i=0; i<_particles.size(); i++) { ControlParticle *cp = &_particles[i]; fprintf(stdout,"A %f,%f \n",cp->density,cp->densityWeight); }
+	avgdw /= (LbmFloat)(_particles.size());
+	//if(motion) { printf("ControlParticle::kernel: avgdw:%f,  kw%f, sp%f \n", avgdw, _kernelWeight, _fluidSpacing); }
+	
+	//if((simtime>=0.) && (simtime != _currTime)) 
+	initTime(simtime, dt);
+
+	if((motion) && (motion->getSize()>0)){
+		ControlParticle *motionp = motion->getParticle(0);
+		//printf("ControlParticle::prepareControl motion: pos[%f,%f,%f] vel[%f,%f,%f] \n", motionp->pos[0], motionp->pos[1], motionp->pos[2], motionp->vel[0], motionp->vel[1], motionp->vel[2] );
+		for(size_t i=0; i<_particles.size(); i++) {
+			ControlParticle *cp = &_particles[i];
+			cp->pos = cp->pos + motionp->pos;
+			cp->vel = cp->vel + motionp->vel;
+			cp->size = cp->size * motionp->size;
+			cp->influence = cp->size * motionp->influence;
+		}
+	}
+	
+	// reset to radiusAtt by default
+	if(_radiusVel==0.) _radiusVel = _radiusAtt;
+	if(_radiusMinMaxd==0.) _radiusMinMaxd = _radiusAtt;
+	if(_radiusMaxd==0.) _radiusMaxd = 2.*_radiusAtt;
+	// has to be radiusVel<radiusAtt<radiusMinMaxd<radiusMaxd
+	if(_radiusVel>_radiusAtt) _radiusVel = _radiusAtt;
+	if(_radiusAtt>_radiusMinMaxd) _radiusAtt = _radiusMinMaxd;
+	if(_radiusMinMaxd>_radiusMaxd) _radiusMinMaxd = _radiusMaxd;
+
+	//printf("ControlParticle::radii vel:%f att:%f min:%f max:%f \n", _radiusVel,_radiusAtt,_radiusMinMaxd,_radiusMaxd);
+	// prepareControl done
+}
+
+void ControlParticles::finishControl(std::vector<ControlForces> &forces, LbmFloat iatt, LbmFloat ivel, LbmFloat imaxd) {
+
+	//const LbmFloat iatt  = this->getInfluenceAttraction() * this->getCurrTimestep();
+	//const LbmFloat ivel  = this->getInfluenceVelocity();
+	//const LbmFloat imaxd = this->getInfluenceMaxdist() * this->getCurrTimestep();
+	// prepare for usage
+	iatt  *= this->getCurrTimestep();
+	ivel  *= 1.; // not necessary!
+	imaxd *= this->getCurrTimestep();
+
+	// skip when size=0
+	for(int i=0; i<(int)forces.size(); i++) {
+		if(DEBUG_MODVEL) fprintf(stdout, "CPFORGF %d , wf:%f,f:%f,%f,%f  ,  v:%f,%f,%f   \n",i, forces[i].weightAtt, forces[i].forceAtt[0],forces[i].forceAtt[1],forces[i].forceAtt[2],   forces[i].forceVel[0], forces[i].forceVel[1], forces[i].forceVel[2] );
+		LbmFloat cfweight = forces[i].weightAtt; // always normalize
+		if((cfweight!=0.)&&(iatt!=0.)) {
+			// multiple kernels, normalize - note this does not normalize in d>r/2 region
+			if(ABS(cfweight)>1.) { cfweight = 1.0/cfweight; }
+			// multiply iatt afterwards to allow stronger force
+			cfweight *= iatt;
+			forces[i].forceAtt *= cfweight;
+		} else {
+			forces[i].weightAtt =  0.;
+			forces[i].forceAtt = LbmVec(0.);
+		}
+
+		if( (cfweight==0.) && (imaxd>0.) && (forces[i].maxDistance>0.) ) {
+			forces[i].forceMaxd *= imaxd;
+		} else {
+			forces[i].maxDistance=  0.;
+			forces[i].forceMaxd = LbmVec(0.);
+		}
+
+		LbmFloat cvweight = forces[i].weightVel; // always normalize
+		if(cvweight>0.) {
+			forces[i].forceVel /= cvweight;
+			forces[i].compAv /= cvweight;
+			// now modify cvweight, and write back
+			// important, cut at 1 - otherwise strong vel. influences...
+			if(cvweight>1.) { cvweight = 1.; }
+			// thus cvweight is in the range of 0..influenceVelocity, currently not normalized by numCParts
+			cvweight *= ivel;
+			if(cvweight<0.) cvweight=0.; if(cvweight>1.) cvweight=1.;
+			// LBM, FIXME todo use relaxation factor
+			//pvel = (cvel*0.5 * cvweight) + (pvel * (1.0-cvweight)); 
+			forces[i].weightVel = cvweight;
+
+			//errMsg("COMPAV","i"<<i<<" compav"<<forces[i].compAv<<" forcevel"<<forces[i].forceVel<<" ");
+		} else {
+			forces[i].weightVel = 0.;
+			if(forces[i].maxDistance==0.) forces[i].forceVel = LbmVec(0.);
+			forces[i].compAvWeight = 0.;
+			forces[i].compAv = LbmVec(0.);
+		}
+		if(DEBUG_MODVEL) fprintf(stdout, "CPFINIF %d , wf:%f,f:%f,%f,%f  ,  v:%f,%f,%f   \n",i, forces[i].weightAtt, forces[i].forceAtt[0],forces[i].forceAtt[1],forces[i].forceAtt[2],  forces[i].forceVel[0],forces[i].forceVel[1],forces[i].forceVel[2] );
+	}
+
+	// unused...
+	if(DEBUG_MODVEL) fprintf(stdout,"MFC iatt:%f,%f ivel:%f,%f ifmd:%f,%f \n", iatt,_radiusAtt, ivel,_radiusVel, imaxd, _radiusMaxd);
+	//for(size_t i=0; i<_particles.size(); i++) { ControlParticle *cp = &_particles[i]; fprintf(stdout," %f,%f,%f ",cp->density,cp->densityWeight, (1.0 - (12.0*cp->densityWeight))); }
+	//fprintf(stdout,"\n\nCP DONE \n\n\n");
+}
+
+
+// --------------------------------------------------------------------------
+// calculate forces at given position, and modify velocity
+// according to timestep
+void ControlParticles::calculateCpInfluenceOpt(ControlParticle *cp, LbmVec fluidpos, LbmVec fluidvel, ControlForces *force, LbmFloat fillFactor) {
+	// dont reset, only add...
+	// test distance, simple squared distance reject
+	const LbmFloat cpfo = _radiusAtt*cp->size;
+
+	LbmVec posDelta;
+	if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f bef fw:%f, f:%f,%f,%f  , vw:%f, v:%f,%f,%f   \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+ 	posDelta	= cp->pos - fluidpos;
+#if LBMDIM==2 && (CP_PROJECT2D==1)
+	posDelta[2] = 0.; // project to xy plane, z-velocity should already be gone...
+#endif
+
+	const LbmFloat distsqr = posDelta[0]*posDelta[0]+posDelta[1]*posDelta[1]+posDelta[2]*posDelta[2];
+	if(DEBUG_MODVEL) fprintf(stdout, " Pd at %f,%f,%f d%f   \n",posDelta[0],posDelta[1],posDelta[2], distsqr);
+	// cut at influence=0.5 , scaling not really makes sense
+	if(cpfo*cpfo < distsqr) {
+		/*if(cp->influence>0.5) {
+			if(force->weightAtt == 0.) {
+				if(force->maxDistance*force->maxDistance > distsqr) {
+				const LbmFloat dis = sqrtf((float)distsqr);
+				const LbmFloat sc = dis-cpfo;
+				force->maxDistance = dis;
+				force->forceMaxd = (posDelta)*(sc/dis);
+				}
+			} } */
+		return;
+	}
+	force->weightAtt += 1e-6; // for distance
+	force->maxDistance = 0.; // necessary for SPH?
+
+	const LbmFloat pdistance = MAGNITUDE(posDelta);
+	LbmFloat pdistinv = 0.;
+	if(ABS(pdistance)>0.) pdistinv = 1./pdistance;
+	posDelta *= pdistinv;
+
+	LbmFloat falloffAtt = 0.; //CPKernel::kernel(cpfo * 1.0, pdistance);
+	const LbmFloat qac = pdistance / cpfo ;
+	if (qac < 1.0){ // return 0.;
+		if(qac < 0.5) falloffAtt =  1.0f;
+		else         falloffAtt = (1.0f - qac) * 2.0f;
+	}
+
+	// vorticity force:
+	// - //LbmVec forceVort; 
+	// - //CROSS(forceVort, posDelta, cp->rotaxis);
+	// - //NORMALIZE(forceVort);
+	// - if(falloffAtt>1.0) falloffAtt=1.0;
+
+#if (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+	// fillFactor *= 2.0 *0.75 * pdistance; // 2d>3d sampling
+#endif // (CP_PROJECT2D==1) && (defined(MAIN_2D) || LBMDIM==2)
+	cp->density += falloffAtt * fillFactor;
+	force->forceAtt += posDelta *cp->densityWeight *cp->influence; 
+	force->weightAtt += falloffAtt*cp->densityWeight *cp->influence;
+	
+	LbmFloat falloffVel = 0.; //CPKernel::kernel(cpfo * 1.0, pdistance);
+	const LbmFloat cpfv = _radiusVel*cp->size;
+	if(cpfv*cpfv < distsqr) { return; }
+	const LbmFloat qvc = pdistance / cpfo ;
+	//if (qvc < 1.0){ 
+		//if(qvc < 0.5) falloffVel =  1.0f;
+		//else         falloffVel = (1.0f - qvc) * 2.0f;
+	//}
+	falloffVel = 1.-qvc;
+
+	LbmFloat pvWeight; // = (1.0-cp->densityWeight) * _currTimestep * falloffVel;
+	pvWeight = falloffVel *cp->influence; // std, without density influence
+	//pvWeight *= (1.0-cp->densityWeight); // use inverse density weight
+	//pvWeight *=      cp->densityWeight; // test, use density weight
+	LbmVec modvel(0.);
+	modvel += cp->vel * pvWeight;
+	//pvWeight = 1.; modvel = partVel; // DEBUG!?
+
+	if(pvWeight>0.) {
+		force->forceVel += modvel;
+		force->weightVel += pvWeight;
+
+		cp->avgVelWeight += falloffVel;
+		cp->avgVel += fluidvel;
+	} 
+	if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f aft fw:%f, f:%f,%f,%f  , vw:%f, v:%f,%f,%f   \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+	return;
+}
+
+void ControlParticles::calculateMaxdForce(ControlParticle *cp, LbmVec fluidpos, ControlForces *force) {
+	if(force->weightAtt != 0.) return; // maxd force off
+	if(cp->influence <= 0.5) return;   // ignore
+
+	LbmVec posDelta;
+	//if(DEBUG_MODVEL) fprintf(stdout, "CP at %f,%f,%f bef fw:%f, f:%f,%f,%f  , vw:%f, v:%f,%f,%f   \n",fluidpos[0],fluidpos[1],fluidpos[2], force->weightAtt, force->forceAtt[0], force->forceAtt[1], force->forceAtt[2], force->weightVel, force->forceVel[0], force->forceVel[1], force->forceVel[2]);
+ 	posDelta	= cp->pos - fluidpos;
+#if LBMDIM==2 && (CP_PROJECT2D==1)
+	posDelta[2] = 0.; // project to xy plane, z-velocity should already be gone...
+#endif
+
+	// dont reset, only add...
+	// test distance, simple squared distance reject
+	const LbmFloat distsqr = posDelta[0]*posDelta[0]+posDelta[1]*posDelta[1]+posDelta[2]*posDelta[2];
+	
+	// closer cp found
+	if(force->maxDistance*force->maxDistance < distsqr) return;
+	
+	const LbmFloat dmin = _radiusMinMaxd*cp->size;
+	if(distsqr<dmin*dmin) return; // inside min
+	const LbmFloat dmax = _radiusMaxd*cp->size;
+	if(distsqr>dmax*dmax) return; // outside
+
+
+	if(DEBUG_MODVEL) fprintf(stdout, " Pd at %f,%f,%f d%f   \n",posDelta[0],posDelta[1],posDelta[2], distsqr);
+	// cut at influence=0.5 , scaling not really makes sense
+	const LbmFloat dis = sqrtf((float)distsqr);
+	//const LbmFloat sc = dis - dmin;
+	const LbmFloat sc = (dis-dmin)/(dmax-dmin); // scale from 0-1
+	force->maxDistance = dis;
+	force->forceMaxd = (posDelta/dis) * sc;
+	//debug errMsg("calculateMaxdForce","pos"<<fluidpos<<" dis"<<dis<<" sc"<<sc<<" dmin"<<dmin<<" maxd"<< force->maxDistance <<" fmd"<<force->forceMaxd );
+	return;
+}
+
diff --git a/intern/elbeem/intern/controlparticles.h b/intern/elbeem/intern/controlparticles.h
new file mode 100644
index 00000000000..5b8683a4f7d
--- /dev/null
+++ b/intern/elbeem/intern/controlparticles.h
@@ -0,0 +1,297 @@
+// --------------------------------------------------------------------------
+//
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the 
+// GNU General Public License. See the file COPYING for details.  
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+// control particle classes
+//
+// --------------------------------------------------------------------------
+
+#ifndef CONTROLPARTICLES_H
+#define CONTROLPARTICLES_H
+
+// indicator for LBM inclusion
+//#ifndef LBMDIM
+
+//#include <NxFoundation.h>
+//#include <vector>
+//class MultisphGUI;
+//#define NORMALIZE(a) a.normalize()
+//#define MAGNITUDE(a) a.magnitude()
+//#define CROSS(a,b,c) a.cross(b,c)
+//#define ABS(a) (a>0. ? (a) : -(a))
+//#include "cpdefines.h"
+
+//#else // LBMDIM
+
+// use compatibility defines
+//#define NORMALIZE(a) normalize(a)
+//#define MAGNITUDE(a) norm(a)
+//#define CROSS(a,b,c) a=cross(b,c)
+
+//#endif // LBMDIM
+
+#define MAGNITUDE(a) norm(a)
+
+// math.h compatibility
+#define CP_PI ((LbmFloat)3.14159265358979323846)
+
+// project 2d test cases onto plane?
+// if not, 3d distance is used for 2d sim as well
+#define CP_PROJECT2D 1
+
+
+// default init for mincpdist, ControlForces::maxDistance
+#define CPF_MAXDINIT 10000.
+
+// storage of influence for a fluid cell/particle in lbm/sph
+class ControlForces
+{
+public:
+	ControlForces() { };
+	~ControlForces() {};
+
+	// attraction force
+	LbmFloat weightAtt;
+	LbmVec forceAtt;
+	// velocity influence
+	LbmFloat weightVel;
+	LbmVec forceVel;
+	// maximal distance influence, 
+	// first is max. distance to first control particle
+	// second attraction strength
+	LbmFloat maxDistance;
+	LbmVec forceMaxd;
+
+	LbmFloat compAvWeight;
+	LbmVec compAv;
+
+	void resetForces() {
+		weightAtt = weightVel = 0.;
+		maxDistance = CPF_MAXDINIT;
+		forceAtt = forceVel = forceMaxd = LbmVec(0.,0.,0.);
+		compAvWeight=0.; compAv=LbmVec(0.);
+	};
+};
+
+
+// single control particle
+class ControlParticle
+{
+public:
+	ControlParticle() { reset(); };
+	~ControlParticle() {};
+
+	// control parameters
+	
+	// position
+	LbmVec pos;
+	// size (influences influence radius)
+	LbmFloat size;
+	// overall strength of influence
+	LbmFloat influence;
+	// rotation axis
+	LbmVec rotaxis;
+
+	// computed values
+
+	// velocity
+	LbmVec vel;
+	// computed density
+	LbmFloat density;
+	LbmFloat densityWeight;
+
+	LbmVec avgVel;
+	LbmVec avgVelAcc;
+	LbmFloat avgVelWeight;
+
+	// init all zero / defaults
+	void reset();
+};
+
+
+// container for a particle configuration at time t
+class ControlParticleSet
+{
+public:
+
+	// time of particle set
+	LbmFloat time;
+	// particle positions
+	std::vector<ControlParticle> particles;
+
+};
+
+
+// container & management of control particles
+class ControlParticles
+{
+public:
+	ControlParticles();
+	~ControlParticles();
+
+	// reset datastructures for next influence step
+	// if motion object is given, particle 1 of second system is used for overall 
+	// position and speed offset
+	void prepareControl(LbmFloat simtime, LbmFloat dt, ControlParticles *motion);
+	// post control operations
+	void finishControl(std::vector<ControlForces> &forces, LbmFloat iatt, LbmFloat ivel, LbmFloat imaxd);
+	// recalculate 
+	void calculateKernelWeight();
+
+	// calculate forces at given position, and modify velocity
+	// according to timestep (from initControl)
+	void calculateCpInfluenceOpt (ControlParticle *cp, LbmVec fluidpos, LbmVec fluidvel, ControlForces *force, LbmFloat fillFactor);
+	void calculateMaxdForce      (ControlParticle *cp, LbmVec fluidpos, ControlForces *force);
+
+	// no. of particles
+	inline int getSize() { return (int)_particles.size(); }
+	int getTotalSize();
+	// get particle [i]
+	inline ControlParticle* getParticle(int i){ return &_particles[i]; }
+
+	// set influence parameters
+	void setInfluenceTangential(LbmFloat set) { _influenceTangential=set; }
+	void setInfluenceAttraction(LbmFloat set) { _influenceAttraction=set; }
+	void setInfluenceMaxdist(LbmFloat set)    { _influenceMaxdist=set; }
+	// calculate for delta t
+	void setInfluenceVelocity(LbmFloat set, LbmFloat dt);
+	// get influence parameters
+	inline LbmFloat getInfluenceAttraction()    { return _influenceAttraction; }
+	inline LbmFloat getInfluenceTangential()    { return _influenceTangential; }
+	inline LbmFloat getInfluenceVelocity()      { return _influenceVelocity; }
+	inline LbmFloat getInfluenceMaxdist()       { return _influenceMaxdist; }
+	inline LbmFloat getCurrTimestep()           { return _currTimestep; }
+
+	void setRadiusAtt(LbmFloat set)       { _radiusAtt=set; }
+	inline LbmFloat getRadiusAtt()        { return _radiusAtt; }
+	void setRadiusVel(LbmFloat set)       { _radiusVel=set; }
+	inline LbmFloat getRadiusVel()        { return _radiusVel; }
+	void setRadiusMaxd(LbmFloat set)      { _radiusMaxd=set; }
+	inline LbmFloat getRadiusMaxd()       { return _radiusMaxd; }
+	void setRadiusMinMaxd(LbmFloat set)   { _radiusMinMaxd=set; }
+	inline LbmFloat getRadiusMinMaxd()    { return _radiusMinMaxd; }
+
+	LbmFloat getControlTimStart();
+	LbmFloat getControlTimEnd();
+
+	// set/get characteristic length (and inverse)
+	void setCharLength(LbmFloat set)      { _charLength=set; _charLengthInv=1./_charLength; }
+	inline LbmFloat getCharLength()       { return _charLength;}
+	inline LbmFloat getCharLengthInv()    { return _charLengthInv;}
+
+	// set init parameters
+	void setInitTimeScale(LbmFloat set)  { _initTimeScale = set; };
+	void setInitMirror(string set)  { _initMirror = set; };
+	string getInitMirror()          { return _initMirror; };
+
+	void setLastOffset(LbmVec set) { _initLastPartOffset = set; };
+	void setLastScale(LbmVec set)  { _initLastPartScale = set; };
+	void setOffset(LbmVec set) { _initPartOffset = set; };
+	void setScale(LbmVec set)  { _initPartScale = set; };
+
+	// set/get cps params
+	void setCPSWith(LbmFloat set)       { mCPSWidth = set; };
+	void setCPSTimestep(LbmFloat set)   { mCPSTimestep = set; };
+	void setCPSTimeStart(LbmFloat set)  { mCPSTimeStart = set; };
+	void setCPSTimeEnd(LbmFloat set)    { mCPSTimeEnd = set; };
+	void setCPSMvmWeightFac(LbmFloat set) { mCPSWeightFac = set; };
+
+	LbmFloat getCPSWith()       { return mCPSWidth; };
+	LbmFloat getCPSTimestep()   { return mCPSTimestep; };
+	LbmFloat getCPSTimeStart()  { return mCPSTimeStart; };
+	LbmFloat getCPSTimeEnd()    { return mCPSTimeEnd; };
+	LbmFloat getCPSMvmWeightFac() { return mCPSWeightFac; };
+
+	void setDebugInit(int set)       { mDebugInit = set; };
+
+	// set init parameters
+	void setFluidSpacing(LbmFloat set)  { _fluidSpacing = set; };
+
+	// load positions & timing from text file
+	int initFromTextFile(string filename);
+	int initFromTextFileOld(string filename);
+	// load positions & timing from gzipped binary file
+	int initFromBinaryFile(string filename);
+	int initFromMVCMesh(string filename);
+	// init an example test case
+	int initExampleSet();
+
+	// init for a given time
+	void initTime(LbmFloat t, LbmFloat dt);
+
+	// blender test init
+	void initBlenderTest();
+
+protected:
+	// sets influence params
+	friend class MultisphGUI;
+
+	// tangential and attraction influence
+	LbmFloat _influenceTangential, _influenceAttraction;
+	// direct velocity influence
+	LbmFloat _influenceVelocity;
+	// maximal distance influence
+	LbmFloat _influenceMaxdist;
+
+	// influence radii
+	LbmFloat _radiusAtt, _radiusVel, _radiusMinMaxd, _radiusMaxd;
+
+	// currently valid time & timestep
+	LbmFloat _currTime, _currTimestep;
+	// all particles
+	std::vector<ControlParticle> _particles;
+
+	// particle sets
+	std::vector<ControlParticleSet> mPartSets;
+
+	// additional parameters for initing particles
+	LbmFloat _initTimeScale;
+	LbmVec _initPartOffset;
+	LbmVec _initPartScale;
+	LbmVec _initLastPartOffset;
+	LbmVec _initLastPartScale;
+	// mirror particles for loading?
+	string _initMirror;
+
+	// row spacing paramter, e.g. use for approximation of kernel area/volume
+	LbmFloat _fluidSpacing;
+	// save current kernel weight
+	LbmFloat _kernelWeight;
+	// charateristic length in world coordinates for normalizatioon of forces
+	LbmFloat _charLength, _charLengthInv;
+
+
+	/*! do ani mesh CPS */
+	void calculateCPS(string filename);
+	//! ani mesh cps params 
+	ntlVec3Gfx mvCPSStart, mvCPSEnd;
+	gfxReal mCPSWidth, mCPSTimestep;
+	gfxReal mCPSTimeStart, mCPSTimeEnd;
+	gfxReal mCPSWeightFac;
+
+	int mDebugInit;
+
+	
+protected:
+	// apply init transformations
+	void applyTrafos();
+
+	// helper function for init -> swap components everywhere
+	void swapCoords(int a,int b);
+	// helper function for init -> mirror time
+	void mirrorTime();
+
+	// helper, init given array
+	void initTimeArray(LbmFloat t, std::vector<ControlParticle> &parts);
+
+	bool checkPointInside(ntlTree *tree, ntlVec3Gfx org, gfxReal &distance);
+};
+
+
+
+#endif
+
diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp
index b2779f51c3b..23a447a3d64 100644
--- a/intern/elbeem/intern/elbeem.cpp
+++ b/intern/elbeem/intern/elbeem.cpp
@@ -30,6 +30,7 @@ int guiRoiMaxLev=6, guiRoiMinLev=0;
 ntlWorld *gpWorld = NULL;
 
 
+
 // API
 
 // reset elbeemSimulationSettings struct with defaults
diff --git a/intern/elbeem/intern/elbeem_control.cpp b/intern/elbeem/intern/elbeem_control.cpp
new file mode 100644
index 00000000000..272cfa6ec9d
--- /dev/null
+++ b/intern/elbeem/intern/elbeem_control.cpp
@@ -0,0 +1,22 @@
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * Control API header
+ */
+
+#include "elbeem.h"
+#include "elbeem_control.h"
+
+// add mesh as fluidsim object
+int elbeemControlAddSet(struct elbeemControl*) {
+	return 0;
+}
+
+int elbeemControlComputeMesh(struct elbeemMesh) {
+	return 0;
+}
+
diff --git a/intern/elbeem/intern/elbeem_control.h b/intern/elbeem/intern/elbeem_control.h
new file mode 100644
index 00000000000..0d30835a224
--- /dev/null
+++ b/intern/elbeem/intern/elbeem_control.h
@@ -0,0 +1,62 @@
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * Control API header
+ */
+#ifndef ELBEEMCONTROL_API_H
+#define ELBEEMCONTROL_API_H
+
+// a single control particle set
+typedef struct elbeemControl {
+	/* influence forces */
+	float influenceAttraction;
+	float *channelInfluenceAttraction;
+	float channelSizeInfluenceAttraction;
+
+	float influenceVelocity;
+	float *channelInfluenceVelocity;
+	float channelSizeInfluenceVelocity;
+
+	float influenceMaxdist;
+	float *channelInfluenceMaxdist;
+	float channelSizeInfluenceMaxdist;
+
+	/* influence force radii */
+	float radiusAttraction;
+	float *channelRadiusAttraction;
+	float channelSizeRadiusAttraction;
+
+	float radiusVelocity;
+	float *channelRadiusVelocity;
+	float channelSizeRadiusVelocity;
+
+	float radiusMindist;
+	float *channelRadiusMindist;
+	float channelSizeRadiusMindist;
+	float radiusMaxdist;
+	float *channelRadiusMaxdist;
+	float channelSizeRadiusMaxdist;
+
+	/* control particle positions/scale */
+	float offset[3];
+	float *channelOffset;
+	float channelSizeOffset;
+	
+	float scale[3];
+	float *channelScale;
+	float channelSizeScale;
+	
+} elbeemControl;
+
+
+// add mesh as fluidsim object
+int elbeemControlAddSet(struct elbeemControl*);
+
+// sample & track mesh control particles, TODO add return type...
+int elbeemControlComputeMesh(struct elbeemMesh);
+
+#endif // ELBEEMCONTROL_API_H
diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp
index 1b0ba13c707..9925565b85d 100644
--- a/intern/elbeem/intern/isosurface.cpp
+++ b/intern/elbeem/intern/isosurface.cpp
@@ -13,10 +13,6 @@
 #include <algorithm>
 #include <stdio.h>
 
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
-
 // just use default rounding for platforms where its not available
 #ifndef round
 #define round(x) (x)
diff --git a/intern/elbeem/intern/mvmcoords.cpp b/intern/elbeem/intern/mvmcoords.cpp
new file mode 100644
index 00000000000..f8ae11ea71d
--- /dev/null
+++ b/intern/elbeem/intern/mvmcoords.cpp
@@ -0,0 +1,191 @@
+/******************************************************************************
+ *
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the 
+// GNU General Public License. See the file COPYING for details.  
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+ *
+ * Mean Value Mesh Coords class
+ *
+ *****************************************************************************/
+
+#include "mvmcoords.h"
+#include <algorithm>
+using std::vector;
+
+void MeanValueMeshCoords::clear() 
+{
+	mVertices.resize(0);
+	mNumVerts = 0;
+}
+
+void MeanValueMeshCoords::calculateMVMCs(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle> &tris, 
+		vector<ntlVec3Gfx> &points, gfxReal numweights)
+{
+	clear();
+	mvmTransferPoint tds;
+	int mem = 0;
+	int i = 0;
+
+	mNumVerts = (int)reference_vertices.size();
+
+	for (vector<ntlVec3Gfx>::iterator iter = points.begin(); iter != points.end(); ++iter, ++i) {
+		if(i%(points.size()/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Computing weights, points: "<<i<<"/"<<points.size(),5 );
+		tds.lastpos = *iter;
+		tds.weights.resize(0);	// clear
+		computeWeights(reference_vertices, tris, tds, numweights);
+		mem += (int)tds.weights.size();
+		mVertices.push_back(tds);
+	}    
+	int mbmem = mem * sizeof(mvmFloat) / (1024*1024);
+	debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"vertices:"<<mNumVerts<<" points:"<<points.size()<<" weights:"<<mem<<", wmem:"<<mbmem<<"MB ",7 );
+}
+
+// from: mean value coordinates for closed triangular meshes
+// attention: fails if a point is exactly (or very close) to a vertex
+void MeanValueMeshCoords::computeWeights(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle>& tris,
+		mvmTransferPoint& tds, gfxReal numweights)
+{
+	const bool mvmFullDebug=false;
+	//const ntlVec3Gfx cEPS = 1.0e-6;
+	const mvmFloat cEPS = 1.0e-14;
+
+	//mvmFloat d[3], s[3], phi[3],c[3];
+	ntlVec3d u[3],c,d,s,phi;
+	int indices[3];
+
+	for (int i = 0; i < (int)reference_vertices.size(); ++i) {
+		tds.weights.push_back(mvmIndexWeight(i, 0.0));
+	}
+
+	// for each triangle
+	//for (vector<ntlTriangle>::iterator iter = tris.begin(); iter != tris.end();) {
+	for(int t=0; t<(int)tris.size(); t++) {
+
+		for (int i = 0; i < 3; ++i) { //, ++iter) {
+			indices[i] = tris[t].getPoints()[i];
+			u[i] = vec2D(reference_vertices[ indices[i] ]-tds.lastpos);
+			d[i] = normalize(u[i]); //.normalize();        
+			//assert(d[i] != 0.);
+			if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","t"<<t<<" i"<<indices[i] //<<" lp"<<tds.lastpos
+					<<" v"<<reference_vertices[indices[i]]<<" u"<<u[i]<<" ");
+			// on vertex!
+			//? if(d[i]<=0.) continue;
+		}
+		//for (int i = 0; i < 3; ++i) { errMsg("III"," "<<i		<<" i"<<indices[i]<<reference_vertices[ indices[i] ] ); }
+
+		// arcsin is not needed, see paper
+		phi[0] = 2.*asin( (mvmFloat)(0.5* norm(u[1]-u[2]) )  );
+		phi[1] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[2]) )  );
+		phi[2] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[1]) )  );
+		mvmFloat h = (phi[0] + phi[1] + phi[2])*0.5;
+		if (M_PI-h < cEPS) {
+			if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","point on triangle");
+			tds.weights.resize(0);
+			tds.weights.push_back( mvmIndexWeight(indices[0], sin(phi[0])*d[1]*d[2]));
+			tds.weights.push_back( mvmIndexWeight(indices[1], sin(phi[1])*d[0]*d[2]));
+			tds.weights.push_back( mvmIndexWeight(indices[2], sin(phi[2])*d[1]*d[0]));
+			break;
+		}
+		mvmFloat sinh = 2.*sin(h);
+		c[0] = (sinh*sin(h-phi[0]))/(sin(phi[1])*sin(phi[2]))-1.;
+		c[1] = (sinh*sin(h-phi[1]))/(sin(phi[0])*sin(phi[2]))-1.;    
+		c[2] = (sinh*sin(h-phi[2]))/(sin(phi[0])*sin(phi[1]))-1.;   
+		if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","c="<<c<<" phi="<<phi<<" d="<<d);
+		//if (c[0] > 1. || c[0] < 0. || c[1] > 1. || c[1] < 0. || c[2] > 1. || c[2] < 0.) continue;
+
+		s[0] = sqrtf((float)(1.-c[0]*c[0]));
+		s[1] = sqrtf((float)(1.-c[1]*c[1]));
+		s[2] = sqrtf((float)(1.-c[2]*c[2]));
+
+		if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","s");
+		if (s[0] <= cEPS || s[1] <= cEPS || s[2] <= cEPS) {
+			//MSG("position lies outside the triangle on the same plane -> ignore it");
+			continue;
+		}
+		const mvmFloat u0x = u[0][0];
+		const mvmFloat u0y = u[0][1];
+		const mvmFloat u0z = u[0][2];
+		const mvmFloat u1x = u[1][0];
+		const mvmFloat u1y = u[1][1];
+		const mvmFloat u1z = u[1][2];
+		const mvmFloat u2x = u[2][0];
+		const mvmFloat u2y = u[2][1];
+		const mvmFloat u2z = u[2][2];
+		mvmFloat det = u0x*u1y*u2z - u0x*u1z*u2y + u0y*u1z*u2x - u0y*u1x*u2z + u0z*u1x*u2y - u0z*u1y*u2x;
+		//assert(det != 0.);
+		if (det < 0.) {
+			s[0] = -s[0];
+			s[1] = -s[1];
+			s[2] = -s[2];
+		}
+
+		tds.weights[indices[0]].weight += (phi[0]-c[1]*phi[2]-c[2]*phi[1])/(d[0]*sin(phi[1])*s[2]);
+		tds.weights[indices[1]].weight += (phi[1]-c[2]*phi[0]-c[0]*phi[2])/(d[1]*sin(phi[2])*s[0]);
+		tds.weights[indices[2]].weight += (phi[2]-c[0]*phi[1]-c[1]*phi[0])/(d[2]*sin(phi[0])*s[1]);
+		if(mvmFullDebug) { errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[0]<<" o"<<tds.weights[indices[0]].weight);
+		errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[1]<<" o"<<tds.weights[indices[1]].weight);
+		errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[2]<<" o"<<tds.weights[indices[2]].weight);
+		errMsg("MeanValueMeshCoords::computeWeights","\n\n\n"); }
+	}
+
+	//sort weights
+	if((numweights>0.)&& (numweights<1.) ) {
+	//if( ((int)tds.weights.size() > maxNumWeights) && (maxNumWeights > 0) ) {
+	  int maxNumWeights = (int)(tds.weights.size()*numweights);
+		if(maxNumWeights<=0) maxNumWeights = 1;
+		std::sort(tds.weights.begin(), tds.weights.end(), std::greater<mvmIndexWeight>());
+		// only use maxNumWeights-th largest weights
+		tds.weights.resize(maxNumWeights);
+	}
+
+	// normalize weights
+	mvmFloat totalWeight = 0.;
+	for (vector<mvmIndexWeight>::const_iterator witer = tds.weights.begin();
+			witer != tds.weights.end(); ++witer) {
+		totalWeight += witer->weight;
+	}
+	mvmFloat invTotalWeight;
+	if (totalWeight == 0.) {
+		if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","totalWeight == 0");
+		invTotalWeight = 0.0;
+	} else {
+		invTotalWeight = 1.0/totalWeight;
+	}
+
+	for (vector<mvmIndexWeight>::iterator viter = tds.weights.begin();
+			viter != tds.weights.end(); ++viter) {
+		viter->weight *= invTotalWeight;  
+		//assert(finite(viter->weight) != 0);
+		if(!finite(viter->weight)) viter->weight=0.;
+	}
+}
+
+void MeanValueMeshCoords::transfer(vector<ntlVec3Gfx> &vertices, vector<ntlVec3Gfx>& displacements) 
+{
+	displacements.resize(0);
+
+	//debMsgStd("MeanValueMeshCoords::transfer",DM_MSG,"vertices:"<<mNumVerts<<" curr_verts:"<<vertices.size()<<" ",7 );
+	if((int)vertices.size() != mNumVerts) {
+		errMsg("MeanValueMeshCoords::transfer","Different no of verts: "<<vertices.size()<<" vs "<<mNumVerts);
+		return;
+	}
+
+	for (vector<mvmTransferPoint>::iterator titer = mVertices.begin(); titer != mVertices.end(); ++titer) {
+		mvmTransferPoint &tds = *titer;
+		ntlVec3Gfx newpos(0.0);
+
+		for (vector<mvmIndexWeight>::iterator witer = tds.weights.begin();
+				witer != tds.weights.end(); ++witer) {
+			newpos += vertices[witer->index] * witer->weight;
+			//errMsg("transfer","np"<<newpos<<" v"<<vertices[witer->index]<<" w"<< witer->weight);
+		}
+
+		displacements.push_back(newpos);
+		//displacements.push_back(newpos - tds.lastpos);
+		//tds.lastpos = newpos;
+	}
+}
+
diff --git a/intern/elbeem/intern/mvmcoords.h b/intern/elbeem/intern/mvmcoords.h
new file mode 100644
index 00000000000..2b63d31da54
--- /dev/null
+++ b/intern/elbeem/intern/mvmcoords.h
@@ -0,0 +1,77 @@
+/******************************************************************************
+ *
+// El'Beem - the visual lattice boltzmann freesurface simulator
+// All code distributed as part of El'Beem is covered by the version 2 of the 
+// GNU General Public License. See the file COPYING for details.  
+//
+// Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
+//
+ *
+ * Mean Value Mesh Coords class
+ *
+ *****************************************************************************/
+
+#ifndef MVMCOORDS_H
+#define MVMCOORDS_H
+
+#include "utilities.h"
+#include "ntl_ray.h"
+#include <vector>
+#define mvmFloat double
+
+// weight and triangle index
+class mvmIndexWeight {
+	public:
+
+		mvmIndexWeight() : weight(0.0) {}
+
+		mvmIndexWeight(int const& i, mvmFloat const& w) :
+			weight(w), index(i) {}
+
+		// for sorting
+		bool operator> (mvmIndexWeight const& w) const { return this->weight > w.weight; } 
+		bool operator< (mvmIndexWeight const& w) const { return this->weight < w.weight; }
+
+		mvmFloat weight;
+		int index;
+};
+
+// transfer point with weights
+class mvmTransferPoint {
+	public:
+		//! position of transfer point
+		ntlVec3Gfx lastpos;
+		//! triangle weights
+		std::vector<mvmIndexWeight> weights;
+};
+
+
+//! compute mvmcs
+class MeanValueMeshCoords {
+
+	public:
+
+    MeanValueMeshCoords() {}
+    ~MeanValueMeshCoords() {
+        clear();
+    }
+
+    void clear();
+
+    void calculateMVMCs(std::vector<ntlVec3Gfx> &reference_vertices, 
+				std::vector<ntlTriangle> &tris, std::vector<ntlVec3Gfx> &points, gfxReal numweights);
+    
+    void transfer(std::vector<ntlVec3Gfx> &vertices, std::vector<ntlVec3Gfx>& displacements);
+
+	protected:
+
+    void computeWeights(std::vector<ntlVec3Gfx> &reference_vertices, 
+				std::vector<ntlTriangle> &tris, mvmTransferPoint& tds, gfxReal numweights);
+
+    std::vector<mvmTransferPoint> mVertices;
+    int mNumVerts;
+
+};
+
+#endif
+
diff --git a/intern/elbeem/intern/solver_adap.cpp b/intern/elbeem/intern/solver_adap.cpp
index 5616d805232..ef516a578bd 100644
--- a/intern/elbeem/intern/solver_adap.cpp
+++ b/intern/elbeem/intern/solver_adap.cpp
@@ -11,9 +11,7 @@
 #include "solver_relax.h"
 #include "particletracer.h"
 
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
+
 
 /*****************************************************************************/
 //! coarse step functions
diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h
index d46f065adfd..310dd50617d 100644
--- a/intern/elbeem/intern/solver_class.h
+++ b/intern/elbeem/intern/solver_class.h
@@ -105,6 +105,10 @@
 #endif
 #endif
 
+#if LBM_INCLUDE_CONTROL==1
+#include "solver_control.h"
+#endif 
+
 #if LBM_INCLUDE_TESTSOLVERS==1
 #include "solver_test.h"
 #endif // LBM_INCLUDE_TESTSOLVERS==1
@@ -497,6 +501,14 @@ class LbmFsgrSolver :
 		LbmFloat& debRAC(LbmFloat* s,int l);
 #		endif // FSGR_STRICT_DEBUG==1
 
+#		if LBM_INCLUDE_CONTROL==1
+		LbmControlData *mpControl;
+
+		void initCpdata();
+		void handleCpdata();
+		void cpDebugDisplay(int dispset); 
+#		endif // LBM_INCLUDE_CONTROL==1
+
 		bool mUseTestdata;
 #		if LBM_INCLUDE_TESTSOLVERS==1
 		// test functions
@@ -506,10 +518,6 @@ class LbmFsgrSolver :
 		void handleTestdata();
 		void set3dHeight(int ,int );
 
-		void initCpdata();
-		void handleCpdata();
-		void cpDebugDisplay(int dispset);
-
 		int mMpNum,mMpIndex;
 		int mOrgSizeX;
 		LbmFloat mOrgStartX;
diff --git a/intern/elbeem/intern/solver_control.cpp b/intern/elbeem/intern/solver_control.cpp
new file mode 100644
index 00000000000..69f1682e253
--- /dev/null
+++ b/intern/elbeem/intern/solver_control.cpp
@@ -0,0 +1,961 @@
+/******************************************************************************
+ *
+ * El'Beem - the visual lattice boltzmann freesurface simulator
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.  
+ *
+ * Copyright 2003-2008 Nils Thuerey 
+ *
+ * control extensions
+ *
+ *****************************************************************************/
+#include "solver_class.h"
+#include "solver_relax.h"
+#include "particletracer.h"
+
+#include "solver_control.h"
+
+#include "controlparticles.h"
+
+
+
+/******************************************************************************
+ * LbmControlData control set
+ *****************************************************************************/
+
+LbmControlSet::LbmControlSet() :
+	mCparts(NULL), mCpmotion(NULL), mContrPartFile(""), mCpmotionFile(""),
+	mcForceAtt(0.), mcForceVel(0.), mcForceMaxd(0.), 
+	mcRadiusAtt(0.), mcRadiusVel(0.), mcRadiusMind(0.), mcRadiusMaxd(0.), 
+	mcCpScale(1.), mcCpOffset(0.)
+{
+}
+LbmControlSet::~LbmControlSet() {
+	if(mCparts) delete mCparts;
+	if(mCpmotion) delete mCpmotion;
+}
+void LbmControlSet::initCparts() {
+	mCparts = new ControlParticles();
+	mCpmotion = new ControlParticles();
+}
+
+
+
+/******************************************************************************
+ * LbmControlData control
+ *****************************************************************************/
+
+LbmControlData::LbmControlData() : 
+	mSetForceStrength(0.),
+	mCons(), mCpUpdateInterval(16), mCpOutfile(""),
+	mCpForces(), mCpKernel(), mMdKernel(),
+	mDiffVelCon(1.),
+	mDebugCpscale(0.), 
+	mDebugVelScale(0.), 
+	mDebugCompavScale(0.), 
+	mDebugAttScale(0.), 
+	mDebugMaxdScale(0.),
+	mDebugAvgVelScale(0.)
+{
+}
+
+LbmControlData::~LbmControlData() 
+{
+}
+
+
+void LbmControlData::parseControldataAttrList(AttributeList *attr) {
+	// controlpart vars
+	mSetForceStrength = attr->readFloat("tforcestrength", mSetForceStrength,"LbmControlData", "mSetForceStrength", false);
+	//errMsg("tforcestrength set to "," "<<mSetForceStrength);
+	mCpUpdateInterval = attr->readInt("controlparticle_updateinterval", mCpUpdateInterval,"LbmControlData","mCpUpdateInterval", false);
+	// tracer output file
+	mCpOutfile = attr->readString("controlparticle_outfile",mCpOutfile,"LbmControlData","mCpOutfile", false);
+	if(getenv("ELBEEM_CPOUTFILE")) {
+		string outfile(getenv("ELBEEM_CPOUTFILE"));
+		mCpOutfile = outfile;
+		debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPOUTFILE to set mCpOutfile to "<<outfile<<","<<mCpOutfile,7);
+	}
+
+	for(int cpii=0; cpii<10; cpii++) {
+		string suffix("");
+		//if(cpii>0)
+		{  suffix = string("0"); suffix[0]+=cpii; }
+		LbmControlSet *cset;
+		cset = new LbmControlSet();
+		cset->initCparts();
+
+		cset->mContrPartFile = attr->readString("controlparticle"+suffix+"_file",cset->mContrPartFile,"LbmControlData","cset->mContrPartFile", false);
+		if((cpii==0) && (getenv("ELBEEM_CPINFILE")) ) {
+			string infile(getenv("ELBEEM_CPINFILE"));
+			cset->mContrPartFile = infile;
+			debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPINFILE to set mContrPartFile to "<<infile<<","<<cset->mContrPartFile,7);
+		}
+
+		LbmFloat cpvort=0.;
+		cset->mcRadiusAtt =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusatt", 0., "LbmControlData","mcRadiusAtt" );
+		cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
+		cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
+		cset->mCparts->setRadiusAtt(cset->mcRadiusAtt.get(0.));
+		cset->mCparts->setRadiusVel(cset->mcRadiusVel.get(0.));
+
+		// WARNING currently only for first set
+		//if(cpii==0) {
+		cset->mcForceAtt  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_attraction", 0. , "LbmControlData","cset->mcForceAtt", false);
+		cset->mcForceVel  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_velocity",   0. , "LbmControlData","mcForceVel", false);
+		cset->mcForceMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_maxdist",    0. , "LbmControlData","mcForceMaxd", false);
+		cset->mCparts->setInfluenceAttraction(cset->mcForceAtt.get(0.) );
+		// warning - stores temprorarily, value converted to dt dep. factor
+		cset->mCparts->setInfluenceVelocity(cset->mcForceVel.get(0.) , 0.01 ); // dummy dt
+		cset->mCparts->setInfluenceMaxdist(cset->mcForceMaxd.get(0.) );
+		cpvort =  attr->readFloat("controlparticle"+suffix+"_vorticity",   cpvort, "LbmControlData","cpvort", false);
+		cset->mCparts->setInfluenceTangential(cpvort);
+			
+		cset->mcRadiusMind =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmin", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMind", false);
+		cset->mcRadiusMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmax", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMaxd", false);
+		cset->mCparts->setRadiusMinMaxd(cset->mcRadiusMind.get(0.));
+		cset->mCparts->setRadiusMaxd(cset->mcRadiusMaxd.get(0.));
+		//}
+
+		// now local...
+		//LbmVec cpOffset(0.), cpScale(1.);
+		LbmFloat cpTimescale = 1.;
+		string cpMirroring("");
+
+		//cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
+		//cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
+		cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
+		cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
+		cpTimescale =  attr->readFloat("controlparticle"+suffix+"_timescale",  cpTimescale, "LbmControlData","cpTimescale", false);
+		cpMirroring =  attr->readString("controlparticle"+suffix+"_mirror",  cpMirroring, "LbmControlData","cpMirroring", false);
+
+		LbmFloat cpsWidth = cset->mCparts->getCPSWith();
+		cpsWidth =  attr->readFloat("controlparticle"+suffix+"_cpswidth",  cpsWidth, "LbmControlData","cpsWidth", false);
+		LbmFloat cpsDt = cset->mCparts->getCPSTimestep();
+		cpsDt =  attr->readFloat("controlparticle"+suffix+"_cpstimestep",  cpsDt, "LbmControlData","cpsDt", false);
+		LbmFloat cpsTstart = cset->mCparts->getCPSTimeStart();
+		cpsTstart =  attr->readFloat("controlparticle"+suffix+"_cpststart",  cpsTstart, "LbmControlData","cpsTstart", false);
+		LbmFloat cpsTend = cset->mCparts->getCPSTimeEnd();
+		cpsTend =  attr->readFloat("controlparticle"+suffix+"_cpstend",  cpsTend, "LbmControlData","cpsTend", false);
+		LbmFloat cpsMvmfac = cset->mCparts->getCPSMvmWeightFac();
+		cpsMvmfac =  attr->readFloat("controlparticle"+suffix+"_cpsmvmfac",  cpsMvmfac, "LbmControlData","cpsMvmfac", false);
+		cset->mCparts->setCPSWith(cpsWidth);
+		cset->mCparts->setCPSTimestep(cpsDt);
+		cset->mCparts->setCPSTimeStart(cpsTstart);
+		cset->mCparts->setCPSTimeEnd(cpsTend);
+		cset->mCparts->setCPSMvmWeightFac(cpsMvmfac);
+
+		cset->mCparts->setOffset( vec2L(cset->mcCpOffset.get(0.)) );
+		cset->mCparts->setScale( vec2L(cset->mcCpScale.get(0.)) );
+		cset->mCparts->setInitTimeScale( cpTimescale );
+		cset->mCparts->setInitMirror( cpMirroring );
+
+		int mDebugInit = 0;
+		mDebugInit = attr->readInt("controlparticle"+suffix+"_debuginit", mDebugInit,"LbmControlData","mDebugInit", false);
+		cset->mCparts->setDebugInit(mDebugInit);
+
+		// motion particle settings
+		LbmVec mcpOffset(0.), mcpScale(1.);
+		LbmFloat mcpTimescale = 1.;
+		string mcpMirroring("");
+
+		cset->mCpmotionFile = attr->readString("cpmotion"+suffix+"_file",cset->mCpmotionFile,"LbmControlData","mCpmotionFile", false);
+		mcpTimescale =  attr->readFloat("cpmotion"+suffix+"_timescale",  mcpTimescale, "LbmControlData","mcpTimescale", false);
+		mcpMirroring =  attr->readString("cpmotion"+suffix+"_mirror",  mcpMirroring, "LbmControlData","mcpMirroring", false);
+		mcpOffset = vec2L( attr->readVec3d("cpmotion"+suffix+"_offset", vec2P(mcpOffset),"LbmControlData","cpOffset", false) );
+		mcpScale =  vec2L( attr->readVec3d("cpmotion"+suffix+"_scale",  vec2P(mcpScale), "LbmControlData","cpScale", false) );
+
+		cset->mCpmotion->setOffset( vec2L(mcpOffset) );
+		cset->mCpmotion->setScale( vec2L(mcpScale) );
+		cset->mCpmotion->setInitTimeScale( mcpTimescale );
+		cset->mCpmotion->setInitMirror( mcpMirroring );
+
+		if(cset->mContrPartFile.length()>1) {
+			errMsg("LbmControlData","Using control particle set "<<cpii<<" file:"<<cset->mContrPartFile<<" cpmfile:"<<cset->mCpmotionFile<<" mirr:'"<<cset->mCpmotion->getInitMirror()<<"' " );
+			mCons.push_back( cset );
+		} else {
+			delete cset;
+		}
+	}
+
+	// debug, testing - make sure theres at least an empty set
+	if(mCons.size()<1) {
+		mCons.push_back( new LbmControlSet() );
+		mCons[0]->initCparts();
+	}
+
+	// take from first set
+	for(int cpii=1; cpii<(int)mCons.size(); cpii++) {
+		mCons[cpii]->mCparts->setRadiusMinMaxd( mCons[0]->mCparts->getRadiusMinMaxd() );
+		mCons[cpii]->mCparts->setRadiusMaxd(    mCons[0]->mCparts->getRadiusMaxd() );
+		mCons[cpii]->mCparts->setInfluenceAttraction( mCons[0]->mCparts->getInfluenceAttraction() );
+		mCons[cpii]->mCparts->setInfluenceTangential( mCons[0]->mCparts->getInfluenceTangential() );
+		mCons[cpii]->mCparts->setInfluenceVelocity( mCons[0]->mCparts->getInfluenceVelocity() , 0.01 ); // dummy dt
+		mCons[cpii]->mCparts->setInfluenceMaxdist( mCons[0]->mCparts->getInfluenceMaxdist() );
+	}
+	
+	// invert for usage in relax macro
+	mDiffVelCon =  1.-attr->readFloat("cpdiffvelcon",  mDiffVelCon, "LbmControlData","mDiffVelCon", false);
+
+	mDebugCpscale     =  attr->readFloat("cpdebug_cpscale",  mDebugCpscale, "LbmControlData","mDebugCpscale", false);
+	mDebugMaxdScale   =  attr->readFloat("cpdebug_maxdscale",  mDebugMaxdScale, "LbmControlData","mDebugMaxdScale", false);
+	mDebugAttScale    =  attr->readFloat("cpdebug_attscale",  mDebugAttScale, "LbmControlData","mDebugAttScale", false);
+	mDebugVelScale    =  attr->readFloat("cpdebug_velscale",  mDebugVelScale, "LbmControlData","mDebugVelScale", false);
+	mDebugCompavScale =  attr->readFloat("cpdebug_compavscale",  mDebugCompavScale, "LbmControlData","mDebugCompavScale", false);
+	mDebugAvgVelScale =  attr->readFloat("cpdebug_avgvelsc",  mDebugAvgVelScale, "LbmControlData","mDebugAvgVelScale", false);
+}
+
+
+void 
+LbmFsgrSolver::initCpdata()
+{
+	// enable for cps via env. vars
+	//if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){ mUseTestdata=1; }
+
+	// NT blender integration manual test setup
+	if(1) {
+		// manually switch on! if this is zero, nothing is done...
+		mpControl->mSetForceStrength = this->mTForceStrength = 1.;
+		mpControl->mCons.clear();
+
+		// add new set
+		LbmControlSet *cset;
+
+		cset = new LbmControlSet();
+		cset->initCparts();
+		// dont load any file
+		cset->mContrPartFile = string("");
+
+		// set radii for attraction & velocity forces
+		// set strength of the forces
+		// don't set directly! but use channels: 
+		// mcForceAtt, mcForceVel, mcForceMaxd, mcRadiusAtt, mcRadiusVel, mcRadiusMind, mcRadiusMaxd etc.
+
+		// wrong: cset->mCparts->setInfluenceAttraction(1.15); cset->mCparts->setRadiusAtt(1.5);
+		// right, e.g., to init some constant values:
+		cset->mcForceAtt = AnimChannel<float>(0.2);
+		cset->mcRadiusAtt = AnimChannel<float>(0.75);
+		cset->mcForceVel = AnimChannel<float>(0.2);
+		cset->mcRadiusVel = AnimChannel<float>(0.75);
+
+		// this value can be left at 0.5:
+		cset->mCparts->setCPSMvmWeightFac(0.5);
+
+		mpControl->mCons.push_back( cset );
+
+		// instead of reading from file (cset->mContrPartFile), manually init some particles
+		mpControl->mCons[0]->mCparts->initBlenderTest();
+
+		// other values that might be interesting to change:
+		//cset->mCparts->setCPSTimestep(0.02);
+		//cset->mCparts->setCPSTimeStart(0.);
+		//cset->mCparts->setCPSTimeEnd(1.); 
+
+		//mpControl->mDiffVelCon = 1.; // more rigid velocity control, 0 (default) allows more turbulence
+	}
+
+	// control particle -------------------------------------------------------------------------------------
+
+	// init cppf stage, use set 0!
+	if(mCppfStage>0) {
+		if(mpControl->mCpOutfile.length()<1) mpControl->mCpOutfile = string("cpout"); // use getOutFilename !?
+		char strbuf[100];
+		const char *cpFormat = "_d%dcppf%d";
+
+		// initial coarse stage, no input
+		if(mCppfStage==1) {
+		 	mpControl->mCons[0]->mContrPartFile = "";
+		} else {
+			// read from prev stage
+			snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage-1);
+			mpControl->mCons[0]->mContrPartFile = mpControl->mCpOutfile;
+			mpControl->mCons[0]->mContrPartFile += strbuf;
+			mpControl->mCons[0]->mContrPartFile += ".cpart2";
+		}
+
+		snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage);
+		mpControl->mCpOutfile += strbuf;
+	} // */
+	
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+		ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+		// now set with real dt
+		cparts->setInfluenceVelocity( mpControl->mCons[cpssi]->mcForceVel.get(0.), mLevel[mMaxRefine].timestep);
+		cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
+		cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
+		errMsg("LbmControlData","CppfStage "<<mCppfStage<<" in:"<<mpControl->mCons[cpssi]->mContrPartFile<<
+				" out:"<<mpControl->mCpOutfile<<" cl:"<< cparts->getCharLength() );
+
+		// control particle test init
+		if(mpControl->mCons[cpssi]->mCpmotionFile.length()>=1) cpmotion->initFromTextFile(mpControl->mCons[cpssi]->mCpmotionFile);
+		// not really necessary...
+		//? cparts->setFluidSpacing( mLevel[mMaxRefine].nodeSize	); // use grid coords!?
+		//? cparts->calculateKernelWeight();
+		//? debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - motion inited: "<<cparts->getSize() ,10);
+
+		// ensure both are on for env. var settings
+		// when no particles, but outfile enabled, initialize
+		const int lev = mMaxRefine;
+		if((mpParticles) && (mpControl->mCpOutfile.length()>=1) && (cpssi==0)) {
+			// check if auto num
+			if( (mpParticles->getNumInitialParticles()<=1) && 
+					(mpParticles->getNumParticles()<=1) ) { // initParticles done afterwards anyway
+				int tracers = 0;
+				const int workSet = mLevel[lev].setCurr;
+				FSGR_FORIJK_BOUNDS(lev) { 
+					if(RFLAG(lev,i,j,k, workSet)&(CFFluid)) tracers++;
+				}
+				if(LBMDIM==3) tracers /= 8;
+				else          tracers /= 4;
+				mpParticles->setNumInitialParticles(tracers);
+				mpParticles->setDumpTextFile(mpControl->mCpOutfile);
+				debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - set tracers #"<<tracers<<", actual #"<<mpParticles->getNumParticles() ,10);
+			}
+			if(mpParticles->getDumpTextInterval()<=0.) {
+				mpParticles->setDumpTextInterval(mLevel[lev].timestep * mLevel[lev].lSizex);
+				debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - dump delta t not set, using dti="<< mpParticles->getDumpTextInterval()<<", sim dt="<<mLevel[lev].timestep, 5 );
+			}
+			mpParticles->setDumpParts(true); // DEBUG? also dump as particle system
+		}
+
+		if(mpControl->mCons[cpssi]->mContrPartFile.length()>=1) cparts->initFromTextFile(mpControl->mCons[cpssi]->mContrPartFile);
+		cparts->setFluidSpacing( mLevel[lev].nodeSize	); // use grid coords!?
+		cparts->calculateKernelWeight();
+		debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles mCons"<<cpssi<<" - inited, parts:"<<cparts->getTotalSize()<<","<<cparts->getSize()<<" dt:"<<mpParam->getTimestep()<<" control time:"<<cparts->getControlTimStart()<<" to "<<cparts->getControlTimEnd() ,10);
+	} // cpssi
+
+	if(getenv("ELBEEM_CPINFILE")) {
+		this->mTForceStrength = 1.0;
+	}
+	this->mTForceStrength = mpControl->mSetForceStrength;
+	if(mpControl->mCpOutfile.length()>=1) mpParticles->setDumpTextFile(mpControl->mCpOutfile);
+
+	// control particle  init end -------------------------------------------------------------------------------------
+
+	// make sure equiv to solver init
+	if(this->mTForceStrength>0.) { \
+		mpControl->mCpForces.resize( mMaxRefine+1 );
+		for(int lev = 0; lev<=mMaxRefine; lev++) {
+			LONGINT rcellSize = (mLevel[lev].lSizex*mLevel[lev].lSizey*mLevel[lev].lSizez);
+			debMsgStd("LbmFsgrSolver::initControl",DM_MSG,"mCpForces init, lev="<<lev<<" rcs:"<<(int)(rcellSize+4)<<","<<(rcellSize*sizeof(ControlForces)/(1024*1024)), 9 );
+			mpControl->mCpForces[lev].resize( (int)(rcellSize+4) );
+			//for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces.push_back( ControlForces() );
+			for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces[lev][i].resetForces();
+		}
+	} // on?
+
+	debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles #mCons "<<mpControl->mCons.size()<<" done", 6);
+}
+
+
+#define CPODEBUG 0
+//define CPINTER ((int)(mpControl->mCpUpdateInterval))
+
+#define KERN(x,y,z)    mpControl->mCpKernel[ (((z)*cpkarWidth + (y))*cpkarWidth + (x)) ]
+#define MDKERN(x,y,z)  mpControl->mMdKernel[ (((z)*mdkarWidth + (y))*mdkarWidth + (x)) ]
+
+#define BOUNDCHECK(x,low,high)  ( ((x)<low) ? low : (((x)>high) ? high : (x) )  )
+#define BOUNDSKIP(x,low,high)  ( ((x)<low) || ((x)>high) )
+
+void 
+LbmFsgrSolver::handleCpdata()
+{
+	myTime_t cpstart = getTime();
+	int cpChecks=0;
+	int cpInfs=0;
+	//debMsgStd("ControlData::handleCpdata",DM_MSG,"called... "<<this->mTForceStrength,1);
+
+	// add cp influence
+  if((true) && (this->mTForceStrength>0.)) {
+		// ok continue...
+	} // on off
+	else {
+		return;
+	}
+	
+	if((mpControl->mCpUpdateInterval<1) || (this->mStepCnt%mpControl->mCpUpdateInterval==0)) {
+		// do full reinit later on... 
+	}
+	else if(this->mStepCnt>mpControl->mCpUpdateInterval) {
+		// only reinit new cells
+		// TODO !? remove loop dependance!?
+#define NOFORCEENTRY(lev, i,j,k) (LBMGET_FORCE(lev, i,j,k).maxDistance==CPF_MAXDINIT) 
+		// interpolate missing
+		for(int lev=0; lev<=mMaxRefine; lev++) {
+		FSGR_FORIJK_BOUNDS(lev) { 
+			if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFFluid|CFInter) )  
+			//if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFInter) )  
+			//if(0)
+			{ // only check new inter? RFLAG?check
+				if(NOFORCEENTRY(lev, i,j,k)) {
+					//errMsg("CP","FE_MISSING at "<<PRINT_IJK<<" f"<<LBMGET_FORCE(lev, i,j,k).weightAtt<<" md"<<LBMGET_FORCE(lev, i,j,k).maxDistance );
+
+					LbmFloat nbs=0.;
+					ControlForces vals;
+					vals.resetForces(); vals.maxDistance = 0.;
+					for(int l=1; l<this->cDirNum; l++) { 
+						int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+						//errMsg("CP","FE_MISSING check "<<PRINT_VEC(ni,nj,nk)<<" f"<<LBMGET_FORCE(lev, ni,nj,nk).weightAtt<<" md"<<LBMGET_FORCE(lev, ni,nj,nk).maxDistance );
+						if(!NOFORCEENTRY(lev, ni,nj,nk)) {
+							//? vals.weightAtt   += LBMGET_FORCE(lev, ni,nj,nk).weightAtt;
+							//? vals.forceAtt    += LBMGET_FORCE(lev, ni,nj,nk).forceAtt;
+							vals.maxDistance += LBMGET_FORCE(lev, ni,nj,nk).maxDistance;
+							vals.forceMaxd   += LBMGET_FORCE(lev, ni,nj,nk).forceMaxd;
+							vals.weightVel   += LBMGET_FORCE(lev, ni,nj,nk).weightVel;
+							vals.forceVel    += LBMGET_FORCE(lev, ni,nj,nk).forceVel; 
+							// ignore att/compAv/avgVel here for now
+							nbs += 1.;
+						}
+					}
+					if(nbs>0.) {
+						nbs = 1./nbs;
+						//? LBMGET_FORCE(lev, i,j,k).weightAtt   =	vals.weightAtt*nbs;
+						//? LBMGET_FORCE(lev, i,j,k).forceAtt    = vals.forceAtt*nbs;
+						LBMGET_FORCE(lev, i,j,k).maxDistance = vals.maxDistance*nbs;
+						LBMGET_FORCE(lev, i,j,k).forceMaxd   =	vals.forceMaxd*nbs;
+						LBMGET_FORCE(lev, i,j,k).weightVel   =	vals.weightVel*nbs;
+						LBMGET_FORCE(lev, i,j,k).forceVel    = vals.forceVel*nbs;
+					}
+						/*ControlForces *ff = &LBMGET_FORCE(lev, i,j,k);  // DEBUG
+						errMsg("CP","FE_MISSING rec at "<<PRINT_IJK // DEBUG
+								<<" w:"<<ff->weightAtt<<" wa:" <<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
+								<<" v:"<<ff->weightVel<<" wv:" <<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
+								<<" v:"<<ff->maxDistance<<" wv:" <<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  ); // DEBUG */
+					// else errMsg("CP","FE_MISSING rec at "<<PRINT_IJK<<" failed!"); // DEBUG
+					
+				}
+			}
+		}} // ijk, lev
+
+		// mStepCnt > mpControl->mCpUpdateInterval
+		return;
+	} else {
+		// nothing to do ...
+		return;
+	}
+
+	// reset
+	for(int lev=0; lev<=mMaxRefine; lev++) {
+		FSGR_FORIJK_BOUNDS(lev) { LBMGET_FORCE(lev,i,j,k).resetForces(); }
+	}
+	// do setup for coarsest level
+	const int coarseLev = 0;
+	const int fineLev = mMaxRefine;
+
+	// init for current time
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+
+		LbmControlSet *cset = mpControl->mCons[cpssi];
+		cparts->setRadiusAtt(cset->mcRadiusAtt.get(mSimulationTime));
+		cparts->setRadiusVel(cset->mcRadiusVel.get(mSimulationTime));
+		cparts->setInfluenceAttraction(cset->mcForceAtt.get(mSimulationTime) );
+		cparts->setInfluenceMaxdist(cset->mcForceMaxd.get(mSimulationTime) );
+		cparts->setRadiusMinMaxd(cset->mcRadiusMind.get(mSimulationTime));
+		cparts->setRadiusMaxd(cset->mcRadiusMaxd.get(mSimulationTime));
+		cparts->calculateKernelWeight(); // always necessary!?
+		cparts->setOffset( vec2L(cset->mcCpOffset.get(mSimulationTime)) );
+		cparts->setScale(  vec2L(cset->mcCpScale.get(mSimulationTime)) );
+
+		cparts->setInfluenceVelocity( cset->mcForceVel.get(mSimulationTime), mLevel[fineLev].timestep );
+		cparts->setLastOffset( vec2L(cset->mcCpOffset.get(mSimulationTime-mLevel[fineLev].timestep)) );
+		cparts->setLastScale(  vec2L(cset->mcCpScale.get(mSimulationTime-mLevel[fineLev].timestep)) );
+	}
+
+	// check actual values
+	LbmFloat iatt  = mpControl->mCons[0]->mCparts->getInfluenceAttraction();
+	LbmFloat ivel  = mpControl->mCons[0]->mCparts->getInfluenceVelocity();
+	LbmFloat imaxd = mpControl->mCons[0]->mCparts->getInfluenceMaxdist();
+	//errMsg("FINCIT","iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
+	for(int cpssi=1; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		LbmFloat iatt2  = mpControl->mCons[cpssi]->mCparts->getInfluenceAttraction();
+		LbmFloat ivel2  = mpControl->mCons[cpssi]->mCparts->getInfluenceVelocity();
+		LbmFloat imaxd2 = mpControl->mCons[cpssi]->mCparts->getInfluenceMaxdist();
+		if(iatt2 >iatt)   iatt = iatt2;
+		if(ivel2 >ivel)   ivel = ivel2;
+		if(imaxd2>imaxd)  imaxd= imaxd2;
+		//errMsg("FINCIT"," "<<cpssi<<" iatt2="<<iatt2<<" ivel2="<<ivel2<<" imaxd2="<<imaxd<<" NEW "<<" iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
+	}
+
+	if(iatt==0. && ivel==0. && imaxd==0.) {
+		debMsgStd("ControlData::initControl",DM_MSG,"Skipped, all zero...",4);
+		return;
+	}
+	//iatt  = mpControl->mCons[1]->mCparts->getInfluenceAttraction(); //ivel  = mpControl->mCons[1]->mCparts->getInfluenceVelocity(); //imaxd = mpControl->mCons[1]->mCparts->getInfluenceMaxdist(); // TTTTTT
+
+	// do control setup
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+		ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+		// TEST!?
+		bool radmod = false;
+		const LbmFloat minRadSize = mLevel[coarseLev].nodeSize * 1.5;
+		if((cparts->getRadiusAtt()>0.) && (cparts->getRadiusAtt()<minRadSize) && (!radmod) ) {
+			LbmFloat radfac = minRadSize / cparts->getRadiusAtt(); radmod=true;
+			debMsgStd("ControlData::initControl",DM_MSG,"Modified radii att, fac="<<radfac, 7);
+			cparts->setRadiusAtt(cparts->getRadiusAtt()*radfac);
+			cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
+			cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+			cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+		} else if((cparts->getRadiusVel()>0.) && (cparts->getRadiusVel()<minRadSize) && (!radmod) ) {
+			LbmFloat radfac = minRadSize / cparts->getRadiusVel();
+			debMsgStd("ControlData::initControl",DM_MSG,"Modified radii vel, fac="<<radfac, 7);
+			cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
+			cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+			cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+		} else if((cparts->getRadiusMaxd()>0.) && (cparts->getRadiusMaxd()<minRadSize) && (!radmod) ) {
+			LbmFloat radfac = minRadSize / cparts->getRadiusMaxd();
+			debMsgStd("ControlData::initControl",DM_MSG,"Modified radii maxd, fac="<<radfac, 7);
+			cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
+			cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
+		}
+		if(radmod) {
+			debMsgStd("ControlData::initControl",DM_MSG,"Modified radii: att="<< 
+				cparts->getRadiusAtt()<<", vel=" << cparts->getRadiusVel()<<", maxd=" <<
+				cparts->getRadiusMaxd()<<", mind=" << cparts->getRadiusMinMaxd() ,5);
+		}
+
+		cpmotion->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), NULL );
+		cparts->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), cpmotion );
+	}
+
+	// do control...
+	for(int lev=0; lev<=mMaxRefine; lev++) {
+		LbmFloat levVolume = 1.;
+		LbmFloat levForceScale = 1.;
+		for(int ll=lev; ll<mMaxRefine; ll++) {
+			if(LBMDIM==3) levVolume *= 8.;
+			else levVolume *= 4.;
+			levForceScale *= 2.;
+		}
+		errMsg("LbmFsgrSolver::handleCpdata","levVolume="<<levVolume<<" levForceScale="<<levForceScale );
+	//todo: scale velocity, att by level timestep!?
+		
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+		// ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+
+		const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
+		LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
+		LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
+		LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
+#if LBMDIM==2
+		gsz = gsx;
+#endif
+		LbmFloat goffx = mvGeoStart[0];
+		LbmFloat goffy = mvGeoStart[1];
+		LbmFloat goffz = mvGeoStart[2];
+
+		//const LbmFloat cpwIncFac = 2.0;
+		// max to two thirds of domain size
+		const int cpw = MIN( mLevel[lev].lSizex/3, MAX( (int)( cparts->getRadiusAtt()  /gsx) +1  , 2) ); // normal kernel, att,vel
+		const int cpkarWidth = 2*cpw+1;
+		mpControl->mCpKernel.resize(cpkarWidth* cpkarWidth* cpkarWidth);
+		ControlParticle cpt; cpt.reset();
+		cpt.pos = LbmVec( (gsx*(LbmFloat)cpw)+goffx, (gsy*(LbmFloat)cpw)+goffy, (gsz*(LbmFloat)cpw)+goffz );  // optimize?
+		cpt.density = 0.5; cpt.densityWeight = 0.5;
+#if LBMDIM==3
+		for(int k= 0; k<cpkarWidth; ++k) {
+#else // LBMDIM==3
+		{ int k = cpw;
+#endif 
+			for(int j= 0; j<cpkarWidth; ++j) 
+				for(int i= 0; i<cpkarWidth; ++i) {
+					KERN(i,j,k).resetForces();
+					//LbmFloat dx = i-cpw; LbmFloat dy = j-cpw; LbmFloat dz = k-cpw;
+					//LbmVec dv = ( LbmVec(dx,dy,dz) );
+					//LbmFloat dl = norm( dv ); //LbmVec dir = dv / dl;
+					LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
+					cparts->calculateCpInfluenceOpt( &cpt, pos, LbmVec(0,0,0), &KERN(i,j,k)  ,1. );
+					/*if((CPODEBUG)&&(k==cpw)) errMsg("kern"," at "<<PRINT_IJK<<" pos"<<pos<<" cpp"<<cpt.pos 
+						<<" wf:"<<KERN(i,j,k).weightAtt<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceAtt[0],KERN(i,j,k).forceAtt[1],KERN(i,j,k).forceAtt[2] )  
+						<<" wf:"<<KERN(i,j,k).weightVel<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceVel[0],KERN(i,j,k).forceVel[1],KERN(i,j,k).forceVel[2] )  
+						<<" wf:"<<KERN(i,j,k).maxDistance<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceMaxd[0],KERN(i,j,k).forceMaxd[1],KERN(i,j,k).forceMaxd[2] )  ); // */
+					KERN(i,j,k).weightAtt *= 2.;
+					KERN(i,j,k).forceAtt *= 2.;
+					//KERN(i,j,k).forceAtt[1] *= 2.; KERN(i,j,k).forceAtt[2] *= 2.;
+					KERN(i,j,k).weightVel *= 2.;
+					KERN(i,j,k).forceVel *= 2.;
+					//KERN(i,j,k).forceVel[1] *= 2.; KERN(i,j,k).forceVel[2] *= 2.;
+				}
+		}
+
+		if(CPODEBUG) errMsg("cpw"," = "<<cpw<<" f"<< cparts->getRadiusAtt()<<" gsx"<<gsx<<" kpw"<<cpkarWidth); // DEBUG
+		// first cp loop - add att and vel forces
+		for(int cppi=0; cppi<cparts->getSize(); cppi++) {
+			ControlParticle *cp = cparts->getParticle(cppi);
+			if(cp->influence<=0.) continue;
+			const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
+			const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
+			int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
+			/*if( ((LBMDIM==3)&&(BOUNDSKIP(cpk - cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
+				((LBMDIM==3)&&(BOUNDSKIP(cpk + cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
+				BOUNDSKIP(cpj - cpwsm, 0, mLevel[lev].lSizey ) ||
+				BOUNDSKIP(cpj + cpwsm, 0, mLevel[lev].lSizey ) ||
+				BOUNDSKIP(cpi - cpwsm, 0, mLevel[lev].lSizex ) ||
+				BOUNDSKIP(cpi + cpwsm, 0, mLevel[lev].lSizex ) ) {
+				continue;
+				} // */
+			int is,ie,js,je,ks,ke;
+			ks = BOUNDCHECK(cpk - cpw, getForZMinBnd(), getForZMaxBnd(lev) );
+			ke = BOUNDCHECK(cpk + cpw, getForZMinBnd(), getForZMaxBnd(lev) );
+			js = BOUNDCHECK(cpj - cpw, 0, mLevel[lev].lSizey );
+			je = BOUNDCHECK(cpj + cpw, 0, mLevel[lev].lSizey );
+			is = BOUNDCHECK(cpi - cpw, 0, mLevel[lev].lSizex );
+			ie = BOUNDCHECK(cpi + cpw, 0, mLevel[lev].lSizex );
+			if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
+			if(CPODEBUG) errMsg("cppft","i"<<cppi<<" cpw"<<cpw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
+			cpInfs++;
+
+			for(int k= ks; k<ke; ++k) {
+				for(int j= js; j<je; ++j) {
+
+					CellFlagType *pflag = &RFLAG(lev,is,j,k, mLevel[lev].setCurr);
+					ControlForces *kk = &KERN( is-cpi+cpw, j-cpj+cpw, k-cpk+cpw);
+					ControlForces *ff = &LBMGET_FORCE(lev,is,j,k); 
+					pflag--; kk--; ff--;
+
+					for(int i= is; i<ie; ++i) {
+						// first cp loop (att,vel)
+						pflag++; kk++; ff++;
+
+						//add weight for bnd cells
+						const LbmFloat pwforce = kk->weightAtt;
+						// control particle mod,
+						// dont add multiple CFFluid fsgr boundaries
+						if(lev==mMaxRefine) {
+							//if( ( ((*pflag)&(CFFluid )) && (lev==mMaxRefine) ) ||
+									//( ((*pflag)&(CFGrNorm)) && (lev <mMaxRefine) ) ) {
+							if((*pflag)&(CFFluid|CFUnused)) {
+								// check not fromcoarse?
+								cp->density += levVolume* kk->weightAtt; // old CFFluid
+							} else if( (*pflag) & (CFEmpty) ) {  
+								cp->density -= levVolume* 0.5; 
+							} else { //if( ((*pflag) & (CFBnd)) ) {  
+								cp->density -= levVolume* 0.2;  // penalty
+							}
+						} else {
+							//if((*pflag)&(CFGrNorm)) {
+								//cp->density += levVolume* kk->weightAtt; // old CFFluid
+							//} 
+						}
+						//else if(!((*pflag) & (CFUnused)) ) {  cp->density -= levVolume* 0.2; } // penalty
+
+						if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
+						{
+
+							cpChecks++;
+							//const LbmFloat pwforce = kk->weightAtt;
+							LbmFloat pwvel = kk->weightVel;
+							if((pwforce==0.)&&(pwvel==0.)) { continue; }
+							ff->weightAtt += 1e-6; // for distance
+
+							if(pwforce>0.) {
+								ff->weightAtt += pwforce *cp->densityWeight *cp->influence;
+								ff->forceAtt += kk->forceAtt *levForceScale *cp->densityWeight *cp->influence;
+
+								// old fill handling here
+								const int workSet =mLevel[lev].setCurr;
+								LbmFloat ux=0., uy=0., uz=0.;
+								FORDF1{  
+									const LbmFloat dfn = QCELL(lev, i,j,k, workSet, l);
+									ux  += (this->dfDvecX[l]*dfn);
+									uy  += (this->dfDvecY[l]*dfn); 
+									uz  += (this->dfDvecZ[l]*dfn); 
+								}
+								// control particle mod
+								cp->avgVelWeight += levVolume*pwforce;
+								cp->avgVelAcc += LbmVec(ux,uy,uz) * levVolume*pwforce;
+							}
+
+							if(pwvel>0.) {
+								// TODO make switch? vel.influence depends on density weight... 
+								// (reduced lowering with 0.75 factor)
+								pwvel *=  cp->influence *(1.-0.75*cp->densityWeight);
+								// control particle mod
+								// todo use Omega instead!?
+								ff->forceVel += cp->vel*levVolume*pwvel * velLatticeScale; // levVolume?
+								ff->weightVel += levVolume*pwvel; // levVolume?
+								ff->compAv += cp->avgVel*levVolume*pwvel; // levVolume?
+								ff->compAvWeight += levVolume*pwvel; // levVolume?
+							}
+
+							if(CPODEBUG) errMsg("cppft","i"<<cppi<<" at "<<PRINT_IJK<<" kern:"<<
+									PRINT_VEC(i-cpi+cpw, j-cpj+cpw, k-cpk+cpw )
+									//<<" w:"<<ff->weightAtt<<" wa:"
+									//<<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
+									//<<" v:"<<ff->weightVel<<" wv:"
+									//<<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
+									//<<" v:"<<ff->maxDistance<<" wv:"
+									//<<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  
+									);
+						} // celltype
+
+					} // ijk
+				} // ijk
+			} // ijk
+		} // cpi, end first cp loop (att,vel)
+		debMsgStd("LbmFsgrSolver::handleCpdata",DM_MSG,"Force cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
+	} //cpssi
+	} // lev
+
+	// second loop
+	for(int lev=0; lev<=mMaxRefine; lev++) {
+		LbmFloat levVolume = 1.;
+		LbmFloat levForceScale = 1.;
+		for(int ll=lev; ll<mMaxRefine; ll++) {
+			if(LBMDIM==3) levVolume *= 8.;
+			else levVolume *= 4.;
+			levForceScale *= 2.;
+		}
+	// prepare maxd forces
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+
+		// WARNING copied from above!
+		const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
+		LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
+		LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
+		LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
+#if LBMDIM==2
+		gsz = gsx;
+#endif
+		LbmFloat goffx = mvGeoStart[0];
+		LbmFloat goffy = mvGeoStart[1];
+		LbmFloat goffz = mvGeoStart[2];
+
+		//const LbmFloat cpwIncFac = 2.0;
+		const int mdw = MIN( mLevel[lev].lSizex/2, MAX( (int)( cparts->getRadiusMaxd() /gsx) +1  , 2) ); // wide kernel, md
+		const int mdkarWidth = 2*mdw+1;
+		mpControl->mMdKernel.resize(mdkarWidth* mdkarWidth* mdkarWidth);
+		ControlParticle cpt; cpt.reset();
+		cpt.density = 0.5; cpt.densityWeight = 0.5;
+		cpt.pos = LbmVec( (gsx*(LbmFloat)mdw)+goffx, (gsy*(LbmFloat)mdw)+goffy, (gsz*(LbmFloat)mdw)+goffz );  // optimize?
+#if LBMDIM==3
+		for(int k= 0; k<mdkarWidth; ++k) {
+#else // LBMDIM==3
+		{ int k = mdw;
+#endif 
+			for(int j= 0; j<mdkarWidth; ++j) 
+				for(int i= 0; i<mdkarWidth; ++i) {
+					MDKERN(i,j,k).resetForces();
+					LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
+					cparts->calculateMaxdForce( &cpt, pos, &MDKERN(i,j,k)  );
+				}
+		}
+
+		// second cpi loop, maxd forces
+		if(cparts->getInfluenceMaxdist()>0.) {
+			for(int cppi=0; cppi<cparts->getSize(); cppi++) {
+				ControlParticle *cp = cparts->getParticle(cppi);
+				if(cp->influence<=0.) continue;
+				const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
+				const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
+				int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
+
+				int is,ie,js,je,ks,ke;
+				ks = BOUNDCHECK(cpk - mdw, getForZMinBnd(), getForZMaxBnd(lev) );
+				ke = BOUNDCHECK(cpk + mdw, getForZMinBnd(), getForZMaxBnd(lev) );
+				js = BOUNDCHECK(cpj - mdw, 0, mLevel[lev].lSizey );
+				je = BOUNDCHECK(cpj + mdw, 0, mLevel[lev].lSizey );
+				is = BOUNDCHECK(cpi - mdw, 0, mLevel[lev].lSizex );
+				ie = BOUNDCHECK(cpi + mdw, 0, mLevel[lev].lSizex );
+				if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
+				if(CPODEBUG) errMsg("cppft","i"<<cppi<<" mdw"<<mdw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
+				cpInfs++;
+
+				for(int k= ks; k<ke; ++k)
+				 for(int j= js; j<je; ++j) {
+					CellFlagType *pflag = &RFLAG(lev,is-1,j,k, mLevel[lev].setCurr);
+					for(int i= is; i<ie; ++i) {
+						// second cpi loop, maxd forces
+						pflag++;
+						if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
+						{
+							cpChecks++;
+							ControlForces *ff = &LBMGET_FORCE(lev,i,j,k); 
+							if(ff->weightAtt == 0.) {
+								ControlForces *kk = &MDKERN( i-cpi+mdw, j-cpj+mdw, k-cpk+mdw);
+								const LbmFloat pmdf = kk->maxDistance;
+								if((ff->maxDistance > pmdf) || (ff->maxDistance<0.))
+									ff->maxDistance = pmdf;
+								ff->forceMaxd = kk->forceMaxd;
+								// todo use Omega instead!?
+								ff->forceVel = cp->vel* velLatticeScale;
+							}
+						} // celltype
+				} } // ijk
+			} // cpi, md loop 
+		} // maxd inf>0 */
+
+
+		debMsgStd("ControlData::initControl",DM_MSG,"Maxd cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
+	} //cpssi
+
+	// normalize, only done once for the whole array
+	mpControl->mCons[0]->mCparts->finishControl( mpControl->mCpForces[lev], iatt,ivel,imaxd );
+
+	} // lev loop
+
+	myTime_t cpend = getTime(); 
+	debMsgStd("ControlData::handleCpdata",DM_MSG,"Time for cpgrid generation:"<< getTimeString(cpend-cpstart)<<", checks:"<<cpChecks<<" infs:"<<cpInfs<<" " ,8);
+
+	// warning, may return  before
+}
+
+#if LBM_USE_GUI==1
+
+#define USE_GLUTILITIES
+#include "../gui/gui_utilities.h"
+
+void LbmFsgrSolver::cpDebugDisplay(int dispset)
+{
+	for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
+		ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
+		//ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
+		// display cp parts
+		const bool cpCubes = false;
+		const bool cpDots = true;
+		const bool cpCpdist = true;
+		const bool cpHideIna = true;
+		glShadeModel(GL_FLAT);
+		glDisable( GL_LIGHTING ); // dont light lines
+
+		// dot influence
+		if((mpControl->mDebugCpscale>0.) && cpDots) {
+			glPointSize(mpControl->mDebugCpscale * 8.);
+			glBegin(GL_POINTS);
+			for(int i=0; i<cparts->getSize(); i++) {
+				if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+				ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
+				//LbmFloat halfsize = 0.5; 
+				LbmFloat scale = cparts->getParticle(i)->densityWeight;
+				//glColor4f( scale,scale,scale,scale );
+				glColor4f( 0.,scale,0.,scale );
+				glVertex3f( org[0],org[1],org[2] );
+				//errMsg("lbmDebugDisplay","CP "<<i<<" at "<<org); // DEBUG
+			}
+			glEnd();
+		}
+
+		// cp positions
+		if((mpControl->mDebugCpscale>0.) && cpDots) {
+			glPointSize(mpControl->mDebugCpscale * 3.);
+			glBegin(GL_POINTS); 
+			glColor3f( 0,1,0 );
+		}
+		for(int i=0; i<cparts->getSize(); i++) {
+			if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+			ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
+			LbmFloat halfsize = 0.5; 
+			LbmFloat scale = cparts->getRadiusAtt() * cparts->getParticle(i)->densityWeight;
+			if(cpCubes){	glLineWidth( 1 ); 
+				glColor3f( 1,1,1 );
+				ntlVec3Gfx s = org-(halfsize * (scale)); 
+				ntlVec3Gfx e = org+(halfsize * (scale)); 
+				drawCubeWire( s,e ); }
+			if((mpControl->mDebugCpscale>0.) && cpDots) {
+				glVertex3f( org[0],org[1],org[2] );
+			}
+		}
+		if(cpDots) glEnd();
+
+		if(mpControl->mDebugAvgVelScale>0.) {
+			const float scale = mpControl->mDebugAvgVelScale;
+
+			glColor3f( 1.0,1.0,1 );
+			glBegin(GL_LINES); 
+			for(int i=0; i<cparts->getSize(); i++) {
+				if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
+				ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
+
+				//errMsg("CPAVGVEL","i"<<i<<" pos"<<org<<" av"<<cparts->getParticle(i)->avgVel);// DEBUG
+				float dx = cparts->getParticle(i)->avgVel[0];
+				float dy = cparts->getParticle(i)->avgVel[1];
+				float dz = cparts->getParticle(i)->avgVel[2];
+				dx *= scale; dy *= scale; dz *= scale;
+				glVertex3f( org[0],org[1],org[2] );
+				glVertex3f( org[0]+dx,org[1]+dy,org[2]+dz );
+			}
+			glEnd();
+		} // */
+
+		if( (LBMDIM==2) && (cpCpdist) ) {
+			
+			// debug, for use of e.g. LBMGET_FORCE LbmControlData *mpControl = this;
+#			define TESTGET_FORCE(lev,i,j,k)   mpControl->mCpForces[lev][ ((k*mLevel[lev].lSizey)+j)*mLevel[lev].lSizex+i ]
+	
+			glBegin(GL_LINES);
+			//const int lev=0; 
+			for(int lev=0; lev<=mMaxRefine; lev++) {
+			FSGR_FORIJK_BOUNDS(lev) {
+					LbmVec pos = LbmVec( 
+						((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex) * ((LbmFloat)i+0.5) + mvGeoStart[0], 
+						((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey) * ((LbmFloat)j+0.5) + mvGeoStart[1], 
+						((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez) * ((LbmFloat)k+0.5) + mvGeoStart[2]  );
+					if(LBMDIM==2) pos[2] = ((mvGeoEnd[2]-mvGeoStart[2])*0.5 + mvGeoStart[2]);
+
+					if((mpControl->mDebugMaxdScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt<=0.) )
+					if(TESTGET_FORCE(lev,i,j,k).maxDistance>=0.) 
+					if(TESTGET_FORCE(lev,i,j,k).maxDistance<CPF_MAXDINIT ) {
+						const float scale = mpControl->mDebugMaxdScale*10001.;
+						float dx = TESTGET_FORCE(lev,i,j,k).forceMaxd[0];
+						float dy = TESTGET_FORCE(lev,i,j,k).forceMaxd[1];
+						float dz = TESTGET_FORCE(lev,i,j,k).forceMaxd[2];
+						dx *= scale; dy *= scale; dz *= scale;
+						glColor3f( 0,1,0 );
+						glVertex3f( pos[0],pos[1],pos[2] );
+						glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+					} // */
+					if((mpControl->mDebugAttScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt>0.)) {
+						const float scale = mpControl->mDebugAttScale*100011.;
+						float dx = TESTGET_FORCE(lev,i,j,k).forceAtt[0];
+						float dy = TESTGET_FORCE(lev,i,j,k).forceAtt[1];
+						float dz = TESTGET_FORCE(lev,i,j,k).forceAtt[2];
+						dx *= scale; dy *= scale; dz *= scale;
+						glColor3f( 1,0,0 );
+						glVertex3f( pos[0],pos[1],pos[2] );
+						glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+					} // */
+					// why check maxDistance?
+					if((mpControl->mDebugVelScale>0.) && (TESTGET_FORCE(lev,i,j,k).maxDistance+TESTGET_FORCE(lev,i,j,k).weightVel>0.)) {
+						float scale = mpControl->mDebugVelScale*1.;
+						float wvscale = TESTGET_FORCE(lev,i,j,k).weightVel;
+						float dx = TESTGET_FORCE(lev,i,j,k).forceVel[0];
+						float dy = TESTGET_FORCE(lev,i,j,k).forceVel[1];
+						float dz = TESTGET_FORCE(lev,i,j,k).forceVel[2];
+						scale *= wvscale;
+						dx *= scale; dy *= scale; dz *= scale;
+						glColor3f( 0.2,0.2,1 );
+						glVertex3f( pos[0],pos[1],pos[2] );
+						glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+					} // */
+					if((mpControl->mDebugCompavScale>0.) && (TESTGET_FORCE(lev,i,j,k).compAvWeight>0.)) {
+						const float scale = mpControl->mDebugCompavScale*1.;
+						float dx = TESTGET_FORCE(lev,i,j,k).compAv[0];
+						float dy = TESTGET_FORCE(lev,i,j,k).compAv[1];
+						float dz = TESTGET_FORCE(lev,i,j,k).compAv[2];
+						dx *= scale; dy *= scale; dz *= scale;
+						glColor3f( 0.2,0.2,1 );
+						glVertex3f( pos[0],pos[1],pos[2] );
+						glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
+					} // */
+			} // att,maxd
+			}
+			glEnd();
+		}
+	} // cpssi
+
+	//fprintf(stderr,"BLA\n");
+	glEnable( GL_LIGHTING ); // dont light lines
+	glShadeModel(GL_SMOOTH);
+}
+
+#else // LBM_USE_GUI==1
+void LbmFsgrSolver::cpDebugDisplay(int dispset) { }
+#endif // LBM_USE_GUI==1
+
+
diff --git a/intern/elbeem/intern/solver_control.h b/intern/elbeem/intern/solver_control.h
new file mode 100644
index 00000000000..ea4ae74e6a6
--- /dev/null
+++ b/intern/elbeem/intern/solver_control.h
@@ -0,0 +1,187 @@
+/******************************************************************************
+ *
+ * El'Beem - the visual lattice boltzmann freesurface simulator
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2006 Nils Thuerey
+ *
+ * testing extensions
+ *
+ *****************************************************************************/
+
+
+#ifndef LBM_TESTCLASS_H
+#define LBM_TESTCLASS_H
+
+//class IsoSurface;
+class ParticleObject;
+class ControlParticles;
+class ControlForces;
+
+//#define NUMGRIDS 2
+//#define MAXNUMSWS 10
+
+// farfield modes
+#define FARF_3DONLY -1
+#define FARF_BOTH    0
+#define FARF_SWEONLY 1
+// dont reuse 3d vars/init
+#define FARF_SEPSWE  2
+
+// relaxation macros for solver_relax.h
+#if LBM_INCLUDE_CONTROL!=1
+
+// defined in relax.h
+
+#else // LBM_INCLUDE_TESTSOLVERS!=1
+
+// WARNING has to match controlparts.h
+#define CPF_ENTRIES     12
+#define CPF_FORCE       0
+#define CPF_VELWEIGHT   3
+#define CPF_VELOCITY    4
+#define CPF_FORCEWEIGHT 7
+#define CPF_MINCPDIST   8
+#define CPF_MINCPDIR    9
+
+#include "controlparticles.h"
+
+// get force entry, set=0 is unused anyway
+#define LBMGET_FORCE(lev, i,j,k)  mpControl->mCpForces[lev][ (LBMGI(lev,i,j,k,0)) ]
+
+// debug mods off...
+// same as in src/solver_relax.h!
+#define __PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+	ux += (grav)[0]; \
+	uy += (grav)[1]; \
+	uz += (grav)[2]; 
+
+//void testMaxdmod(int i, int j,int k, LbmFloat &ux,LbmFloat &uy,LbmFloat &uz,ControlForces &ff);
+#if LBMDIM==3
+#define MAXDGRAV \
+			if(myforce->forceMaxd[0]*ux+myforce->forceMaxd[1]*uy<LBM_EPSILON) { \
+				ux = v2w*myforce->forceVel[0]+ v2wi*ux;  \
+				uy = v2w*myforce->forceVel[1]+ v2wi*uy; }  \
+			/* movement inverse to g? */ \
+			if((uz>LBM_EPSILON)&&(uz>myforce->forceVel[2])) { \
+				uz = v2w*myforce->forceVel[2]+ v2wi*uz; } 
+#else // LBMDIM==3
+#define MAXDGRAV \
+			if(myforce->forceMaxd[0]*ux<LBM_EPSILON) { \
+				ux = v2w*myforce->forceVel[0]+ v2wi*ux; } \
+			/* movement inverse to g? */ \
+			if((uy>LBM_EPSILON)&&(uy>myforce->forceVel[1])) { \
+				uy = v2w*myforce->forceVel[1]+ v2wi*uy; } 
+#endif // LBMDIM==3
+
+// debug modifications of collide vars (testing)
+// requires: lev,i,j,k
+#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+	LbmFloat attforce = 1.; \
+	if(this->mTForceStrength>0.) { \
+		ControlForces* myforce = &LBMGET_FORCE(lev,i,j,k); \
+		const LbmFloat vf = myforce->weightAtt;\
+		const LbmFloat vw = myforce->weightVel;\
+		if(vf!=0.) { attforce = MAX(0., 1.-vf);  /* TODO FIXME? use ABS(vf) for repulsion force? */ \
+			ux += myforce->forceAtt[0]; \
+			uy += myforce->forceAtt[1]; \
+			uz += myforce->forceAtt[2]; \
+			\
+		} else if(( myforce->maxDistance>0.) && ( myforce->maxDistance<CPF_MAXDINIT)) {\
+			const LbmFloat v2w = mpControl->mCons[0]->mCparts->getInfluenceMaxdist() * \
+				(myforce->maxDistance-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()) / (mpControl->mCons[0]->mCparts->getRadiusMaxd()-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()); \
+			const LbmFloat v2wi = 1.-v2w; \
+			if(v2w>0.){ MAXDGRAV; \
+				/* errMsg("ERRMDTT","at "<<PRINT_IJK<<" maxd="<<myforce->maxDistance<<", newu"<<PRINT_VEC(ux,uy,uz)<<", org"<<PRINT_VEC(oux,ouy,ouz)<<", fv"<<myforce->forceVel<<" " );  */ \
+			}\
+		} \
+		if(vw>0.) { \
+			const LbmFloat vwi = 1.-vw;\
+			const LbmFloat vwd = mpControl->mDiffVelCon;\
+			ux += vw*(myforce->forceVel[0]-myforce->compAv[0] + vwd*(myforce->compAv[0]-ux) ); \
+			uy += vw*(myforce->forceVel[1]-myforce->compAv[1] + vwd*(myforce->compAv[1]-uy) ); \
+			uz += vw*(myforce->forceVel[2]-myforce->compAv[2] + vwd*(myforce->compAv[2]-uz) ); \
+			/*  TODO test!? modify smooth vel by influence of force for each lbm step, to account for force update only each N steps */ \
+			myforce->compAv = (myforce->forceVel*vw+ myforce->compAv*vwi); \
+		} \
+	} \
+	ux += (grav)[0]*attforce; \
+	uy += (grav)[1]*attforce; \
+	uz += (grav)[2]*attforce; \
+ 	/* end PRECOLLIDE_MODS */
+
+#define TEST_IF_CHECK \
+		if((!iffilled)&&(LBMGET_FORCE(lev,i,j,k).weightAtt!=0.)) { \
+			errMsg("TESTIFFILL"," at "<<PRINT_IJK<<" "<<mass<<" "<<rho); \
+			iffilled = true; \
+			if(mass<rho*1.0) mass = rho*1.0; myfrac = 1.0; \
+		}
+
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
+
+
+// a single set of control particles and params
+class LbmControlSet {
+	public:
+		LbmControlSet();
+		~LbmControlSet();
+		void initCparts();
+
+		// control particles
+		ControlParticles *mCparts; 
+		// control particle overall motion (for easier manual generation)
+		ControlParticles *mCpmotion;
+		// cp data file
+		string mContrPartFile;
+		string mCpmotionFile;
+		// cp debug displau
+		LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
+
+		// params
+		AnimChannel<float> mcForceAtt;
+		AnimChannel<float> mcForceVel;
+		AnimChannel<float> mcForceMaxd;
+
+		AnimChannel<float> mcRadiusAtt;
+		AnimChannel<float> mcRadiusVel;
+		AnimChannel<float> mcRadiusMind;
+		AnimChannel<float> mcRadiusMaxd;
+
+		AnimChannel<ntlVec3f> mcCpScale;
+		AnimChannel<ntlVec3f> mcCpOffset;
+};
+		
+
+
+// main control data storage
+class LbmControlData 
+{
+	public:
+		LbmControlData();
+		virtual ~LbmControlData();
+
+		// control data
+
+		// contorl params
+		void parseControldataAttrList(AttributeList *attr);
+
+		// control strength, set for solver interface
+		LbmFloat mSetForceStrength;
+		// cp vars
+		std::vector<LbmControlSet*> mCons;
+		// update interval
+		int mCpUpdateInterval;
+		// output
+		string mCpOutfile;
+		// control particle precomputed influence
+		std::vector< std::vector<ControlForces> > mCpForces;
+		std::vector<ControlForces> mCpKernel;
+		std::vector<ControlForces> mMdKernel;
+		// activate differential velcon
+		LbmFloat mDiffVelCon;
+
+		// cp debug displau
+		LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
+};
+
+#endif // LBM_TESTCLASS_H
diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp
index c953d2f47da..6dd654f57ff 100644
--- a/intern/elbeem/intern/solver_init.cpp
+++ b/intern/elbeem/intern/solver_init.cpp
@@ -328,7 +328,10 @@ LbmFsgrSolver::LbmFsgrSolver() :
 	mInit2dYZ(false),
 	mForceTadapRefine(-1), mCutoff(-1)
 {
-  // not much to do here... 
+#if LBM_INCLUDE_CONTROL==1
+	mpControl = new LbmControlData();
+#endif
+
 #if LBM_INCLUDE_TESTSOLVERS==1
 	mpTest = new LbmTestdata();
 	mMpNum = mMpIndex = 0;
@@ -488,6 +491,10 @@ void LbmFsgrSolver::parseAttrList()
 	mSimulationTime += starttimeskip;
 	if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
 
+#if LBM_INCLUDE_CONTROL==1
+	mpControl->parseControldataAttrList(mpSifAttrs);
+#endif
+
 #if LBM_INCLUDE_TESTSOLVERS==1
 	mUseTestdata = 0;
 	mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
@@ -703,7 +710,7 @@ bool LbmFsgrSolver::initializeSolverMemory()
 			memBlockAllocProblem = true;
 		}
 #endif // Mac
-		if(sizeof(void *)==4 && memEstFine>maxDefaultMemChunk) {
+		if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
 			// max memory chunk for 32bit systems 2gig
 			memBlockAllocProblem = true;
 		}
@@ -1264,8 +1271,11 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
 	mInitDone = 1;
 
-#if LBM_INCLUDE_TESTSOLVERS==1
+#if LBM_INCLUDE_CONTROL==1
 	initCpdata();
+#endif // LBM_INCLUDE_CONTROL==1
+
+#if LBM_INCLUDE_TESTSOLVERS==1
 	initTestdata();
 #endif // ELBEEM_PLUGIN!=1
 	// not inited? dont use...
diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp
index 13ebf91b696..616ac0e2a23 100644
--- a/intern/elbeem/intern/solver_main.cpp
+++ b/intern/elbeem/intern/solver_main.cpp
@@ -13,11 +13,6 @@
 #include "loop_tools.h"
 #include <stdlib.h>
 
-#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
-#include <ieeefp.h>
-#endif
-
-
 /*****************************************************************************/
 /*! perform a single LBM step */
 /*****************************************************************************/
@@ -58,7 +53,7 @@ void LbmFsgrSolver::stepMain() {
 
 	// init moving bc's, can change mMaxVlen
 	initMovingObstacles(false);
-#if LBM_INCLUDE_TESTSOLVERS==1
+#if LBM_INCLUDE_CONTROL==1
 	handleCpdata();
 #endif
 
diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h
index 287b73c77b2..db6e537aff3 100644
--- a/intern/elbeem/intern/solver_relax.h
+++ b/intern/elbeem/intern/solver_relax.h
@@ -15,7 +15,8 @@
 #define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
 #endif // FSGR_STRICT_DEBUG==1
 
-#if LBM_INCLUDE_TESTSOLVERS!=1
+// #if LBM_INCLUDE_TESTSOLVERS!=1
+#if LBM_INCLUDE_CONTROL!=1
 
 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
 	ux += (grav)[0]; \
@@ -24,42 +25,9 @@
 
 #define TEST_IF_CHECK 
 
-#else // LBM_INCLUDE_TESTSOLVERS!=1
-// defined in test.h
-
-#define NEWDIRVELMOTEST 0
-#if NEWDIRVELMOTEST==1
-// off for non testing
-#undef PRECOLLIDE_MODS
-#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
-	ux += (grav)[0]; \
-	uy += (grav)[1]; \
-	uz += (grav)[2];  \
-  { \
-		int lev = mMaxRefine, nomb=0; \
-		LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
-		for(int l=1; l<this->cDfNum; l++) {  \
-			if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
-				if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
-					nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
-					nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
-					bcnt += 1.; \
-				}	else { \
-					nomb++; \
-				} \
-			}  \
-		} \
-		if((bcnt>0.)&&(nomb==0)) { \
-			ux = nux/bcnt; \
-			uy = nuy/bcnt; \
-			uz = nuz/bcnt;  \
-		} \
-	} 
-#else // NEWDIRVELMOTEST==1
-// off for non testing
-#endif // NEWDIRVELMOTEST==1
-
-#endif // LBM_INCLUDE_TESTSOLVERS!=1
+#else // LBM_INCLUDE_CONTROL!=1
+// defined in solver_control.h 
+#endif // LBM_INCLUDE_CONTROL!=1
 
 	
 /******************************************************************************