/** \file elbeem/intern/solver_interface.cpp
 *  \ingroup elbeem
 */
/******************************************************************************
 *
 * 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
 *
 * Combined 2D/3D Lattice Boltzmann Interface Class
 * contains stuff to be statically compiled
 * 
 *****************************************************************************/

/* LBM Files */ 
#include "ntl_matrices.h"
#include "solver_interface.h" 
#include "ntl_ray.h"
#include "ntl_world.h"
#include "elbeem.h"




/******************************************************************************
 * Interface Constructor
 *****************************************************************************/
LbmSolverInterface::LbmSolverInterface() :
	mPanic( false ),
  mSizex(10), mSizey(10), mSizez(10), 
  mAllfluid(false), mStepCnt( 0 ),
	mFixMass( 0.0 ),
  mOmega( 1.0 ),
  mGravity(0.0),
	mSurfaceTension( 0.0 ), 
	mBoundaryEast(  (CellFlagType)(CFBnd) ),mBoundaryWest( (CellFlagType)(CFBnd) ),mBoundaryNorth(  (CellFlagType)(CFBnd) ),
	mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop(  (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
  mInitDone( false ),
  mInitDensityGradient( false ),
	mpSifAttrs( NULL ), mpSifSwsAttrs(NULL), mpParam( NULL ), mpParticles(NULL),
	mNumParticlesLost(0), 
	mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
	mDebugVelScale( 0.01 ), mNodeInfoString("+"),
	mvGeoStart(-1.0), mvGeoEnd(1.0), mpSimTrafo(NULL),
	mAccurateGeoinit(0),
	mName("lbm_default") ,
	mpIso( NULL ), mIsoValue(0.499),
	mSilent(false) , 
	mLbmInitId(1) ,
	mpGiTree( NULL ),
	mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ),
	mRefinementDesired(0),
	mOutputSurfacePreview(0), mPreviewFactor(0.25),
	mSmoothSurface(1.0), mSmoothNormals(1.0),
	mIsoSubdivs(1), mPartGenProb(0.),
	mDumpVelocities(false),
	mMarkedCells(), mMarkedCellIndex(0),
	mDomainBound("noslip"), mDomainPartSlipValue(0.1),
	mFarFieldSize(0.), 
	mPartDropMassSub(0.1),  // good default
	mPartUsePhysModel(false),
	mTForceStrength(0.0), 
	mCppfStage(0),
	mDumpRawText(false),
	mDumpRawBinary(false),
	mDumpRawBinaryZip(true)
#if PARALLEL==1
	, mNumOMPThreads(1)
#endif  // PARALLEL==1
{
#if ELBEEM_PLUGIN==1
	if(gDebugLevel<=1) setSilent(true);
#endif
	mpSimTrafo = new ntlMat4Gfx(0.0); 
	mpSimTrafo->initId();

	if(getenv("ELBEEM_RAWDEBUGDUMP")) { 
		debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2);
		mDumpRawText = true;
	}

	if(getenv("ELBEEM_BINDEBUGDUMP")) { 
		debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2);
		mDumpRawBinary = true;
	}
}

LbmSolverInterface::~LbmSolverInterface() 
{ 
	if(mpSimTrafo) delete mpSimTrafo;
}


/******************************************************************************
 * initialize correct grid sizes given a geometric bounding box
 * and desired grid resolutions, all params except maxrefine
 * will be modified
 *****************************************************************************/
void initGridSizes(int &sizex, int &sizey, int &sizez,
		ntlVec3Gfx &geoStart, ntlVec3Gfx &geoEnd, 
		int mMaxRefine, bool parallel) 
{
	// fix size inits to force cubic cells and mult4 level dimensions
	const int debugGridsizeInit = 1;
  if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" "<<geoStart<<","<<geoEnd ,10);

	int maxGridSize = sizex; // get max size
	if(sizey>maxGridSize) maxGridSize = sizey;
	if(sizez>maxGridSize) maxGridSize = sizez;
	LbmFloat maxGeoSize = (geoEnd[0]-geoStart[0]); // get max size
	if((geoEnd[1]-geoStart[1])>maxGeoSize) maxGeoSize = (geoEnd[1]-geoStart[1]);
	if((geoEnd[2]-geoStart[2])>maxGeoSize) maxGeoSize = (geoEnd[2]-geoStart[2]);
	// FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though
	LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize);
  if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Start:"<<geoStart<<" End:"<<geoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10);
	// force grid sizes according to geom. size, rounded
	sizex = (int) ((geoEnd[0]-geoStart[0]) / cellSize +0.5);
	sizey = (int) ((geoEnd[1]-geoStart[1]) / cellSize +0.5);
	sizez = (int) ((geoEnd[2]-geoStart[2]) / cellSize +0.5);
	// match refinement sizes, round downwards to multiple of 4
	int sizeMask = 0;
	int maskBits = mMaxRefine;
	if(parallel==1) maskBits+=2;
	for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }

	// at least size 4 on coarsest level
	int minSize = 4<<mMaxRefine; //(maskBits+2);
	if(sizex<minSize) sizex = minSize;
	if(sizey<minSize) sizey = minSize;
	if(sizez<minSize) sizez = minSize;
	
	sizeMask = ~sizeMask;
  if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask)<<", maskBits:"<<maskBits<<",minSize:"<<minSize ,10);
	sizex &= sizeMask;
	sizey &= sizeMask;
	sizez &= sizeMask;
	// debug
	int nextinc = (~sizeMask)+1;
  debMsgStd("initGridSizes",DM_MSG,"MPTeffMLtest, curr size:"<<PRINT_VEC(sizex,sizey,sizez)<<", next size:"<<PRINT_VEC(sizex+nextinc,sizey+nextinc,sizez+nextinc) ,10);

	// force geom size to match rounded/modified grid sizes
	geoEnd[0] = geoStart[0] + cellSize*(LbmFloat)sizex;
	geoEnd[1] = geoStart[1] + cellSize*(LbmFloat)sizey;
	geoEnd[2] = geoStart[2] + cellSize*(LbmFloat)sizez;
}

void calculateMemreqEstimate( int resx,int resy,int resz, 
		int refine, float farfield,
		double *reqret, double *reqretFine, string *reqstr) {
	// debug estimation?
	const bool debugMemEst = true;
	// COMPRESSGRIDS define is not available here, make sure it matches
	const bool useGridComp = true;
	// make sure we can handle bid numbers here... all double
	double memCnt = 0.0;
	double ddTotalNum = (double)dTotalNum;
	if(reqretFine) *reqretFine = -1.;

	double currResx = (double)resx;
	double currResy = (double)resy;
	double currResz = (double)resz;
	double rcellSize = ((currResx*currResy*currResz) *ddTotalNum);
	memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
	if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"res:"<<PRINT_VEC(currResx,currResy,currResz)<<" rcellSize:"<<rcellSize<<" mc:"<<memCnt, 10);
  if(!useGridComp) {
		memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
		if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
		if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," no-comp, mc:"<<memCnt, 10);
	} else {
		double compressOffset = (double)(currResx*currResy*ddTotalNum*2.0);
		memCnt += (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0));
		if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0));
		if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," w-comp, mc:"<<memCnt, 10);
	}
	for(int i=refine-1; i>=0; i--) {
		currResx /= 2.0;
		currResy /= 2.0;
		currResz /= 2.0;
		rcellSize = ((currResz*currResy*currResx) *ddTotalNum);
		memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
		memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
		if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"refine "<<i<<", mc:"<<memCnt, 10);
	}

	// isosurface memory, use orig res values
	if(farfield>0.) {
		memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
	} else {
		// ignore 3 int slices...
		memCnt += (double)( (              sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
	}
	if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"iso, mc:"<<memCnt, 10);

	// cpdata init check missing...

	double memd = memCnt;
	const char *sizeStr = "";
	const double sfac = 1024.0;
	if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
	if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
	if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
	if(memd>sfac){ memd /= sfac; sizeStr="TB"; }

	// return values
	std::ostringstream ret;
	if(memCnt< 1024.0*1024.0) {
		// show full MBs
		ret << (ceil(memd));
	} else {
		// two digits for anything larger than MB
		ret << (ceil(memd*100.0)/100.0);
	}
	ret	<< " "<< sizeStr;
	*reqret = memCnt;
	*reqstr = ret.str();
	//debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4);
}

void LbmSolverInterface::initDomainTrafo(float *mat) { 
	mpSimTrafo->initArrayCheck(mat); 
}

/*******************************************************************************/
/*! parse a boundary flag string */
CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) {
	string val  = mpSifAttrs->readString(name, "", source, target, needed);
	if(!strcmp(val.c_str(),"")) {
		// no value given...
		return CFEmpty;
	}
	if(!strcmp(val.c_str(),"bnd_no")) {
		return (CellFlagType)( CFBnd );
	}
	if(!strcmp(val.c_str(),"bnd_free")) {
		return (CellFlagType)( CFBnd );
	}
	if(!strcmp(val.c_str(),"fluid")) {
		/* might be used for some in/out flow cases */
		return (CellFlagType)( CFFluid );
	}
	errMsg("LbmSolverInterface::readBoundaryFlagInt","Invalid value '"<<val<<"' " );
# if FSGR_STRICT_DEBUG==1
	errFatal("readBoundaryFlagInt","Strict abort..."<<val, SIMWORLD_INITERROR);
# endif
	return defaultValue;
}

/*******************************************************************************/
/*! parse standard attributes */
void LbmSolverInterface::parseStdAttrList() {
	if(!mpSifAttrs) {
		errFatal("LbmSolverInterface::parseAttrList","mpSifAttrs pointer not initialized!",SIMWORLD_INITERROR);
		return; }

	// st currently unused
	mBoundaryEast  = readBoundaryFlagInt("boundary_east",  mBoundaryEast, "LbmSolverInterface", "mBoundaryEast", false);
	mBoundaryWest  = readBoundaryFlagInt("boundary_west",  mBoundaryWest, "LbmSolverInterface", "mBoundaryWest", false);
	mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmSolverInterface", "mBoundaryNorth", false);
	mBoundarySouth = readBoundaryFlagInt("boundary_south", mBoundarySouth,"LbmSolverInterface", "mBoundarySouth", false);
	mBoundaryTop   = readBoundaryFlagInt("boundary_top",   mBoundaryTop,"LbmSolverInterface", "mBoundaryTop", false);
	mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmSolverInterface", "mBoundaryBottom", false);

	mpSifAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo ); 

	LbmVec sizeVec(mSizex,mSizey,mSizez);
	sizeVec = vec2L( mpSifAttrs->readVec3d("size",  vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) );
	mSizex = (int)sizeVec[0]; 
	mSizey = (int)sizeVec[1]; 
	mSizez = (int)sizeVec[2];
	// param needs size in any case
	// test solvers might not have mpPara, though
	if(mpParam) mpParam->setSize(mSizex, mSizey, mSizez ); 

	mInitDensityGradient = mpSifAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false);
	mIsoValue = mpSifAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false );

	mDebugVelScale = mpSifAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false);
	mNodeInfoString = mpSifAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false );

	mDumpVelocities = mpSifAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false );
	if(getenv("ELBEEM_DUMPVELOCITIES")) {
		int get = atoi(getenv("ELBEEM_DUMPVELOCITIES"));
		if((get==0)||(get==1)) {
			mDumpVelocities = get;
			debMsgStd("LbmSolverInterface::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_DUMPVELOCITIES to set mDumpVelocities to "<<get<<","<<mDumpVelocities,8);
		}
	}
	
	// new test vars
	// mTForceStrength set from test solver 
	mFarFieldSize = mpSifAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false);
	// old compat
	float sizeScale = mpSifAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false);
	if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); }

	mPartDropMassSub = mpSifAttrs->readFloat("part_dropmass", mPartDropMassSub,"LbmSolverInterface", "mPartDropMassSub", false);

	mCppfStage = mpSifAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false);
	mPartGenProb = mpSifAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false);
	mPartUsePhysModel = mpSifAttrs->readBool("partusephysmodel", mPartUsePhysModel,"LbmFsgrSolver", "mPartUsePhysModel", false);
	mIsoSubdivs = mpSifAttrs->readInt("isosubdivs", mIsoSubdivs,"LbmFsgrSolver", "mIsoSubdivs", false);
}


/*******************************************************************************/
/*! geometry initialization */
/*******************************************************************************/

/*****************************************************************************/
/*! init tree for certain geometry init */
void LbmSolverInterface::initGeoTree() {
	if(mpGlob == NULL) { errFatal("LbmSolverInterface::initGeoTree","Requires globals!",SIMWORLD_INITERROR); return; }
	ntlScene *scene = mpGlob->getSimScene();
	mpGiObjects = scene->getObjects();
	mGiObjInside.resize( mpGiObjects->size() );
	mGiObjDistance.resize( mpGiObjects->size() );
	mGiObjSecondDist.resize( mpGiObjects->size() );
	for(size_t i=0; i<mpGiObjects->size(); i++) { 
		if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
		//debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 );
	}
	debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9)
	debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Objects: "<<mpGiObjects->size() ,10);

	if(mpGiTree != NULL) delete mpGiTree;
	char treeFlag = (1<<(this->mLbmInitId+4));
	mpGiTree = new ntlTree( 
# if FSGR_STRICT_DEBUG!=1
			15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
# else // FSGR_STRICT_DEBUG==1
			10, 20,  // reduced/slower debugging values
# endif // FSGR_STRICT_DEBUG==1
			scene, treeFlag );
}

/*****************************************************************************/
/*! destroy tree etc. when geometry init done */
void LbmSolverInterface::freeGeoTree() {
	if(mpGiTree != NULL) {
		delete mpGiTree;
		mpGiTree = NULL;
	}
}


int globGeoInitDebug = 0;
int globGICPIProblems = 0;
/*****************************************************************************/
/*! check for a certain flag type at position org */
bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance, int shootDir) {
	// shift ve ctors to avoid rounding errors
	org += ntlVec3Gfx(0.0001);
	OId = -1;

	// select shooting direction
	ntlVec3Gfx dir = ntlVec3Gfx(1., 0., 0.);
	if(shootDir==1) dir = ntlVec3Gfx(0., 1., 0.);
	else if(shootDir==2) dir = ntlVec3Gfx(0., 0., 1.);
	else if(shootDir==-2) dir = ntlVec3Gfx(0., 0., -1.);
	else if(shootDir==-1) dir = ntlVec3Gfx(0., -1., 0.);
	
	ntlRay ray(org, dir, 0, 1.0, mpGlob);
	bool done = false;
	bool inside = false;

	vector<int> giObjFirstHistSide;
	giObjFirstHistSide.resize( mpGiObjects->size() );
	for(size_t i=0; i<mGiObjInside.size(); i++) { 
		mGiObjInside[i] = 0; 
		mGiObjDistance[i] = -1.0; 
		mGiObjSecondDist[i] = -1.0; 
		giObjFirstHistSide[i] = 0; 
	}
	// if not inside, return distance to first hit
	gfxReal firstHit=-1.0;
	int     firstOId = -1;
	if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org<<" f"<<flags<<" acc"<<mAccurateGeoinit);
	
	if(mAccurateGeoinit) {
		while(!done) {
			// find first inside intersection
			ntlTriangle *triIns = NULL;
			distance = -1.0;
			ntlVec3Gfx normal(0.0);
			if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
			else            mpGiTree->intersect(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[OId]++;
					if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = 1;
					if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
				} else {
					// inside hit
					mGiObjInside[OId]++;
					if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
					if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
					if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = -1;
				}
				norg += normal * getVecEpsilon();
				ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
				// 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;
		for(size_t i=0; i<mGiObjInside.size(); i++) {
			if(mGiObjInside[i]>0) {
				bool mess = false;
				if((mGiObjInside[i]%2)==1) {
					if(giObjFirstHistSide[i] != -1) mess=true;
				} else {
					if(giObjFirstHistSide[i] !=  1) mess=true;
				}
				if(mess) {
					//errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
					globGICPIProblems++;
					mGiObjInside[i]++; // believe first hit side...
				}
			}
		}
		for(size_t i=0; i<mGiObjInside.size(); i++) {
			if(globGeoInitDebug) errMsg("CHIII","i"<<i<<" ins="<<mGiObjInside[i]<<" t"<<mGiObjDistance[i]<<" d"<<distance);
			if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
				if(  (distance<0.0)                             || // first intersection -> good
					  ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
						) {						
					distance = mGiObjDistance[i];
					OId = i;
					inside = true;
				} 
			}
		}
		if(!inside) {
			distance = firstHit;
			OId = firstOId;
		}
		if(globGeoInitDebug) errMsg("CHIII","i"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);

		return inside;
	} else {

		// find first inside intersection
		ntlTriangle *triIns = NULL;
		distance = -1.0;
		ntlVec3Gfx normal(0.0);
		if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
		else            mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
		if(triIns) {
			// check outside intersect
			LbmFloat orientation = dot(normal, dir);
			if(orientation<=0.0) return false;

			OId = triIns->getObjectId();
			return true;
		}
		return false;
	}
}
bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse) {
	// shift ve ctors to avoid rounding errors
	org += ntlVec3Gfx(0.0001); //?
	OId = -1;
	ntlRay ray(org, dir, 0, 1.0, mpGlob);
	//int insCnt = 0;
	bool done = false;
	bool inside = false;
	for(size_t i=0; i<mGiObjInside.size(); i++) { 
		mGiObjInside[i] = 0; 
		mGiObjDistance[i] = -1.0; 
		mGiObjSecondDist[i] = -1.0; 
	}
	// if not inside, return distance to first hit
	gfxReal firstHit=-1.0;
	int     firstOId = -1;
	thinHit = false;
	
	if(mAccurateGeoinit) {
		while(!done) {
			// find first inside intersection
			ntlTriangle *triIns = NULL;
			distance = -1.0;
			ntlVec3Gfx normal(0.0);
			mpGiTree->intersect(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;
					//mGiObjDistance[OId] = -1.0;
					//errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
				} else {
					// inside hit
					//if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
					//errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
					if(mGiObjInside[OId]==1) {
						// second inside hit
						if(mGiObjSecondDist[OId]<0.0) mGiObjSecondDist[OId] = distance;
					}
				}
				mGiObjInside[OId]++;
				// always store first hit for thin obj init
				if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;

				norg += normal * getVecEpsilon();
				ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
				// 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;
				//if(insCnt%2) inside=true;
			}
		}

		distance = -1.0;
		// standard inside check
		for(size_t i=0; i<mGiObjInside.size(); i++) {
			if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
				if(  (distance<0.0)                             || // first intersection -> good
					  ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
						) {						
					distance = mGiObjDistance[i];
					OId = i;
					inside = true;
				} 
			}
		}
		// now check for thin hits
		if(!inside) {
			distance = -1.0;
			for(size_t i=0; i<mGiObjInside.size(); i++) {
				if((mGiObjInside[i]>=2)&&(mGiObjDistance[i]>0.0)&&(mGiObjSecondDist[i]>0.0)&&
					 (mGiObjDistance[i]<1.0*halfCellsize)&&(mGiObjSecondDist[i]<2.0*halfCellsize) ) {
					if(  (distance<0.0)                             || // first intersection -> good
							((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
							) {						
						distance = mGiObjDistance[i];
						OId = i;
						inside = true;
						thinHit = true;
					} 
				}
			}
		}
		if(!inside) {
			// check for hit in this cell, opposite to current dir (only recurse once)
			if(recurse) {
				gfxReal r_distance;
				int r_OId;
				bool ret = geoInitCheckPointInside(org, dir*-1.0, flags, r_OId, r_distance, halfCellsize, thinHit, false);
				if((ret)&&(thinHit)) {
					OId = r_OId;
					distance = 0.0; 
					return true;
				}
			}
		}
		// really no hit...
		if(!inside) {
			distance = firstHit;
			OId = firstOId;
			/*if((mGiObjDistance[OId]>0.0)&&(mGiObjSecondDist[OId]>0.0)) {
				const gfxReal thisdist = mGiObjSecondDist[OId]-mGiObjDistance[OId];
				// dont walk over this cell...
				if(thisdist<halfCellsize) distance-=2.0*halfCellsize;
			} // ? */
		}
		//errMsg("CHIII","i"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);

		return inside;
	} else {

		// find first inside intersection
		ntlTriangle *triIns = NULL;
		distance = -1.0;
		ntlVec3Gfx normal(0.0);
		mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
		if(triIns) {
			// check outside intersect
			LbmFloat orientation = dot(normal, dir);
			if(orientation<=0.0) return false;

			OId = triIns->getObjectId();
			return true;
		}
		return false;
	}
}


/*****************************************************************************/
ntlVec3Gfx LbmSolverInterface::getGeoMaxMovementVelocity(LbmFloat simtime, LbmFloat stepsize) {
	ntlVec3Gfx max(0.0);
	if(mpGlob == NULL) return max;
	// mpGiObjects has to be inited here...
	
	for(int i=0; i< (int)mpGiObjects->size(); i++) {
		//errMsg("LbmSolverInterface::getGeoMaxMovementVelocity","i="<<i<<" "<< (*mpGiObjects)[i]->getName() ); // DEBUG
		if( (*mpGiObjects)[i]->getGeoInitType() & (FGI_FLUID|FGI_MBNDINFLOW) ){
			//ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
			ntlVec3Gfx orgvel = (*mpGiObjects)[i]->calculateMaxVel( simtime, simtime+stepsize );
			if( normNoSqrt(orgvel) > normNoSqrt(max) ) { max = orgvel; } 
			//errMsg("MVT","i"<<i<<" a"<< (*mpGiObjects)[i]->getInitialVelocity(simtime)<<" o"<<orgvel ); // DEBUG
			// TODO check if inflow and simtime 
			ntlVec3Gfx inivel = (*mpGiObjects)[i]->getInitialVelocity(simtime);
			if( normNoSqrt(inivel) > normNoSqrt(max) ) { max = inivel; } 
		}
	}
	errMsg("LbmSolverInterface::getGeoMaxMovementVelocity", "max="<< max ); // DEBUG
	return max;
}




/*******************************************************************************/
/*! cell iteration functions */
/*******************************************************************************/




/*****************************************************************************/
//! add cell to mMarkedCells list
void 
LbmSolverInterface::addCellToMarkedList( CellIdentifierInterface *cid ) {
	for(size_t i=0; i<mMarkedCells.size(); i++) {
		// check if cids alreay in
		if( mMarkedCells[i]->equal(cid) ) return;
		//mMarkedCells[i]->setEnd(false);
	}
	mMarkedCells.push_back( cid );
	//cid->setEnd(true);
}

/*****************************************************************************/
//! marked cell iteration methods
CellIdentifierInterface* 
LbmSolverInterface::markedGetFirstCell( ) {
	if(mMarkedCells.size() > 0){ return mMarkedCells[0]; }
	return NULL;
}

CellIdentifierInterface* 
LbmSolverInterface::markedAdvanceCell() {
	mMarkedCellIndex++;
	if(mMarkedCellIndex>=(int)mMarkedCells.size()) return NULL;
	return mMarkedCells[mMarkedCellIndex];
}

void LbmSolverInterface::markedClearList() {
	// FIXME free cids?
	mMarkedCells.clear();
}

#if PARALLEL==1
void LbmSolverInterface::setNumOMPThreads(int num_threads) {
	mNumOMPThreads = num_threads;
}
#endif  // PARALLEL==1

/*******************************************************************************/
/*! string helper functions */
/*******************************************************************************/



// 32k
string convertSingleFlag2String(CellFlagType cflag) {
	CellFlagType flag = cflag;
	if(flag == CFUnused         ) return string("cCFUnused");
	if(flag == CFEmpty          ) return string("cCFEmpty");      
	if(flag == CFBnd            ) return string("cCFBnd");        
	if(flag == CFBndNoslip      ) return string("cCFBndNoSlip");        
	if(flag == CFBndFreeslip    ) return string("cCFBndFreeSlip");        
	if(flag == CFBndPartslip    ) return string("cCFBndPartSlip");        
	if(flag == CFNoInterpolSrc  ) return string("cCFNoInterpolSrc");
	if(flag == CFFluid          ) return string("cCFFluid");      
	if(flag == CFInter          ) return string("cCFInter");      
	if(flag == CFNoNbFluid      ) return string("cCFNoNbFluid");  
	if(flag == CFNoNbEmpty      ) return string("cCFNoNbEmpty");  
	if(flag == CFNoDelete       ) return string("cCFNoDelete");   
	if(flag == CFNoBndFluid     ) return string("cCFNoBndFluid"); 
	if(flag == CFGrNorm         ) return string("cCFGrNorm");     
	if(flag == CFGrFromFine     ) return string("cCFGrFromFine"); 
	if(flag == CFGrFromCoarse   ) return string("cCFGrFromCoarse");
	if(flag == CFGrCoarseInited ) return string("cCFGrCoarseInited");
	if(flag == CFMbndInflow )     return string("cCFMbndInflow");
	if(flag == CFMbndOutflow )    return string("cCFMbndOutflow");
	if(flag == CFInvalid        ) return string("cfINVALID");

	std::ostringstream mult;
	int val = 0;
	if(flag != 0) {
		while(! (flag&1) ) {
			flag = flag>>1;
			val++;
		}
	} else {
		val = -1;
	}
	if(val>=24) {
		mult << "cfOID_" << (flag>>24) <<"_TYPE";
	} else {
		mult << "cfUNKNOWN_" << val <<"_TYPE";
	}
	return mult.str();
}
	
//! helper function to convert flag to string (for debuggin)
string convertCellFlagType2String( CellFlagType cflag ) {
	int flag = cflag;

	const int jmax = sizeof(CellFlagType)*8;
	bool somefound = false;
	std::ostringstream mult;
	mult << "[";
	for(int j=0; j<jmax ; j++) {
		if(flag& (1<<j)) {
			if(somefound) mult << "|";
			mult << j<<"<"<< convertSingleFlag2String( (CellFlagType)(1<<j) ); // this call should always be _non_-recursive
			somefound = true;
		}
	};
	mult << "]";

	// return concatenated string
	if(somefound) return mult.str();

	// empty?
	return string("[emptyCFT]");
}