- another minor solver update to fix
[blender.git] / intern / elbeem / intern / solver_interface.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * All code distributed as part of El'Beem is covered by the version 2 of the 
5  * GNU General Public License. See the file COPYING for details.
6  * Copyright 2003-2005 Nils Thuerey
7  *
8  * Combined 2D/3D Lattice Boltzmann Interface Class
9  * contains stuff to be statically compiled
10  * 
11  *****************************************************************************/
12
13 /* LBM Files */ 
14 #include "ntl_matrices.h"
15 #include "solver_interface.h" 
16 #include "ntl_ray.h"
17 #include "ntl_world.h"
18 #include "elbeem.h"
19
20
21
22
23 /******************************************************************************
24  * Interface Constructor
25  *****************************************************************************/
26 LbmSolverInterface::LbmSolverInterface() :
27         mPanic( false ),
28   mSizex(10), mSizey(10), mSizez(10), 
29   mAllfluid(false), mStepCnt( 0 ),
30         mFixMass( 0.0 ),
31   mOmega( 1.0 ),
32   mGravity(0.0),
33         mSurfaceTension( 0.0 ), 
34         mBoundaryEast(  (CellFlagType)(CFBnd) ),mBoundaryWest( (CellFlagType)(CFBnd) ),mBoundaryNorth(  (CellFlagType)(CFBnd) ),
35         mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop(  (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
36   mInitDone( false ),
37   mInitDensityGradient( false ),
38         mpAttrs( NULL ), mpParam( NULL ), mpParticles(NULL),
39         mNumParticlesLost(0), 
40         mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
41         mDebugVelScale( 0.01 ), mNodeInfoString("+"),
42         mvGeoStart(-1.0), mvGeoEnd(1.0), mpSimTrafo(NULL),
43         mAccurateGeoinit(0),
44         mName("lbm_default") ,
45         mpIso( NULL ), mIsoValue(0.499),
46         mSilent(false) , 
47         mLbmInitId(1) ,
48         mpGiTree( NULL ),
49         mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ),
50         mRefinementDesired(0),
51         mOutputSurfacePreview(0), mPreviewFactor(0.25),
52         mSmoothSurface(1.0), mSmoothNormals(1.0),
53         mPartGenProb(0.),
54         mDumpVelocities(false),
55         mMarkedCells(), mMarkedCellIndex(0),
56         mDomainBound("noslip"), mDomainPartSlipValue(0.1),
57         mTForceStrength(0.0), mFarFieldSize(0.), mCppfStage(0)
58 {
59 #if ELBEEM_PLUGIN==1
60         if(gDebugLevel<=1) mSilent = true;
61 #endif
62         mpSimTrafo = new ntlMat4Gfx(0.0); 
63         mpSimTrafo->initId();
64 }
65
66 LbmSolverInterface::~LbmSolverInterface() 
67
68         if(mpSimTrafo) delete mpSimTrafo;
69 }
70
71
72 /******************************************************************************
73  * initialize correct grid sizes given a geometric bounding box
74  * and desired grid resolutions, all params except maxrefine
75  * will be modified
76  *****************************************************************************/
77 void initGridSizes(int &sizex, int &sizey, int &sizez,
78                 ntlVec3Gfx &geoStart, ntlVec3Gfx &geoEnd, 
79                 int mMaxRefine, bool parallel) 
80 {
81         // fix size inits to force cubic cells and mult4 level dimensions
82         const int debugGridsizeInit = 0;
83   if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" " ,10);
84
85         int maxGridSize = sizex; // get max size
86         if(sizey>maxGridSize) maxGridSize = sizey;
87         if(sizez>maxGridSize) maxGridSize = sizez;
88         LbmFloat maxGeoSize = (geoEnd[0]-geoStart[0]); // get max size
89         if((geoEnd[1]-geoStart[1])>maxGeoSize) maxGeoSize = (geoEnd[1]-geoStart[1]);
90         if((geoEnd[2]-geoStart[2])>maxGeoSize) maxGeoSize = (geoEnd[2]-geoStart[2]);
91         // FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though
92         LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize);
93   if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Start:"<<geoStart<<" End:"<<geoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10);
94         // force grid sizes according to geom. size, rounded
95         sizex = (int) ((geoEnd[0]-geoStart[0]) / cellSize +0.5);
96         sizey = (int) ((geoEnd[1]-geoStart[1]) / cellSize +0.5);
97         sizez = (int) ((geoEnd[2]-geoStart[2]) / cellSize +0.5);
98         // match refinement sizes, round downwards to multiple of 4
99         int sizeMask = 0;
100         int maskBits = mMaxRefine;
101         if(parallel==1) maskBits+=2;
102         for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }
103
104         // at least size 4 on coarsest level
105         int minSize = 2<<(maskBits+2);
106         if(sizex<minSize) sizex = minSize;
107         if(sizey<minSize) sizey = minSize;
108         if(sizez<minSize) sizez = minSize;
109         
110         sizeMask = ~sizeMask;
111   if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask) ,10);
112         sizex &= sizeMask;
113         sizey &= sizeMask;
114         sizez &= sizeMask;
115
116         // force geom size to match rounded/modified grid sizes
117         geoEnd[0] = geoStart[0] + cellSize*(LbmFloat)sizex;
118         geoEnd[1] = geoStart[1] + cellSize*(LbmFloat)sizey;
119         geoEnd[2] = geoStart[2] + cellSize*(LbmFloat)sizez;
120 }
121
122 void calculateMemreqEstimate( int resx,int resy,int resz, int refine,
123                 double *reqret, string *reqstr) {
124         // debug estimation?
125         const bool debugMemEst = false;
126         // COMPRESSGRIDS define is not available here, make sure it matches
127         const bool useGridComp = true;
128         // make sure we can handle bid numbers here... all double
129         double memCnt = 0.0;
130         double ddTotalNum = (double)dTotalNum;
131
132         double currResx = (double)resx;
133         double currResy = (double)resy;
134         double currResz = (double)resz;
135         double rcellSize = ((currResx*currResy*currResz) *ddTotalNum);
136         memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
137         if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"res:"<<PRINT_VEC(currResx,currResy,currResz)<<" rcellSize:"<<rcellSize<<" mc:"<<memCnt, 10);
138   if(!useGridComp) {
139                 memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
140                 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," no-comp, mc:"<<memCnt, 10);
141         } else {
142                 double compressOffset = (double)(currResx*currResy*ddTotalNum*2.0);
143                 memCnt += (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0));
144                 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," w-comp, mc:"<<memCnt, 10);
145         }
146         for(int i=refine-1; i>=0; i--) {
147                 currResx /= 2.0;
148                 currResy /= 2.0;
149                 currResz /= 2.0;
150                 rcellSize = ((currResz*currResy*currResx) *ddTotalNum);
151                 memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0);
152                 memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0);
153                 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"refine "<<i<<", mc:"<<memCnt, 10);
154         }
155
156         // isosurface memory
157         memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
158         if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"iso, mc:"<<memCnt, 10);
159
160         double memd = memCnt;
161         char *sizeStr = "";
162         const double sfac = 1000.0;
163         if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
164         if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
165         if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
166         if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
167
168         // return values
169         std::ostringstream ret;
170         if(memCnt< 1024.0*1024.0) {
171                 // show full MBs
172                 ret << (ceil(memd));
173         } else {
174                 // two digits for anything larger than MB
175                 ret << (ceil(memd*100.0)/100.0);
176         }
177         ret     << " "<< sizeStr;
178         *reqret = memCnt;
179         *reqstr = ret.str();
180         //debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4);
181 }
182
183 void LbmSolverInterface::initDomainTrafo(float *mat) { 
184         mpSimTrafo->initArrayCheck(mat); 
185 }
186
187 /*******************************************************************************/
188 /*! parse a boundary flag string */
189 CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) {
190         string val  = mpAttrs->readString(name, "", source, target, needed);
191         if(!strcmp(val.c_str(),"")) {
192                 // no value given...
193                 return CFEmpty;
194         }
195         if(!strcmp(val.c_str(),"bnd_no")) {
196                 return (CellFlagType)( CFBnd );
197         }
198         if(!strcmp(val.c_str(),"bnd_free")) {
199                 return (CellFlagType)( CFBnd );
200         }
201         if(!strcmp(val.c_str(),"fluid")) {
202                 /* might be used for some in/out flow cases */
203                 return (CellFlagType)( CFFluid );
204         }
205         errMsg("LbmSolverInterface::readBoundaryFlagInt","Invalid value '"<<val<<"' " );
206 # if FSGR_STRICT_DEBUG==1
207         errFatal("readBoundaryFlagInt","Strict abort..."<<val, SIMWORLD_INITERROR);
208 # endif
209         return defaultValue;
210 }
211
212 /*******************************************************************************/
213 /*! parse standard attributes */
214 void LbmSolverInterface::parseStdAttrList() {
215         if(!mpAttrs) {
216                 errFatal("LbmSolverInterface::parseAttrList","mpAttrs pointer not initialized!",SIMWORLD_INITERROR);
217                 return; }
218
219         // st currently unused
220         //mSurfaceTension  = mpAttrs->readFloat("d_surfacetension",  mSurfaceTension, "LbmSolverInterface", "mSurfaceTension", false);
221         mBoundaryEast  = readBoundaryFlagInt("boundary_east",  mBoundaryEast, "LbmSolverInterface", "mBoundaryEast", false);
222         mBoundaryWest  = readBoundaryFlagInt("boundary_west",  mBoundaryWest, "LbmSolverInterface", "mBoundaryWest", false);
223         mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmSolverInterface", "mBoundaryNorth", false);
224         mBoundarySouth = readBoundaryFlagInt("boundary_south", mBoundarySouth,"LbmSolverInterface", "mBoundarySouth", false);
225         mBoundaryTop   = readBoundaryFlagInt("boundary_top",   mBoundaryTop,"LbmSolverInterface", "mBoundaryTop", false);
226         mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmSolverInterface", "mBoundaryBottom", false);
227
228         mpAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo ); 
229
230         LbmVec sizeVec(mSizex,mSizey,mSizez);
231         sizeVec = vec2L( mpAttrs->readVec3d("size",  vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) );
232         mSizex = (int)sizeVec[0]; 
233         mSizey = (int)sizeVec[1]; 
234         mSizez = (int)sizeVec[2];
235         // param needs size in any case
236         // test solvers might not have mpPara, though
237         if(mpParam) mpParam->setSize(mSizex, mSizey, mSizez ); 
238
239         mInitDensityGradient = mpAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false);
240         mIsoValue = mpAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false );
241
242         mDebugVelScale = mpAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false);
243         mNodeInfoString = mpAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false );
244
245         mDumpVelocities = mpAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false );
246         if(getenv("ELBEEM_DUMPVELOCITIES")) {
247                 int get = atoi(getenv("ELBEEM_DUMPVELOCITIES"));
248                 if((get==0)||(get==1)) {
249                         mDumpVelocities = get;
250                         debMsgStd("LbmSolverInterface::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_DUMPVELOCITIES to set mDumpVelocities to "<<get<<","<<mDumpVelocities,8);
251                 }
252         }
253         
254         // new test vars
255         mTForceStrength = 0.; // set from test solver mpAttrs->readFloat("tforcestrength", mTForceStrength,"LbmSolverInterface", "mTForceStrength", false);
256         mFarFieldSize = mpAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false);
257         // old compat
258         float sizeScale = mpAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false);
259         if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); }
260
261         mCppfStage = mpAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false);
262         mPartGenProb = mpAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false);
263 }
264
265
266 /*******************************************************************************/
267 /*! geometry initialization */
268 /*******************************************************************************/
269
270 /*****************************************************************************/
271 /*! init tree for certain geometry init */
272 void LbmSolverInterface::initGeoTree() {
273         if(mpGlob == NULL) { errFatal("LbmSolverInterface::initGeoTree","Requires globals!",SIMWORLD_INITERROR); return; }
274         ntlScene *scene = mpGlob->getSimScene();
275         mpGiObjects = scene->getObjects();
276         mGiObjInside.resize( mpGiObjects->size() );
277         mGiObjDistance.resize( mpGiObjects->size() );
278         mGiObjSecondDist.resize( mpGiObjects->size() );
279         for(size_t i=0; i<mpGiObjects->size(); i++) { 
280                 if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
281                 debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 );
282         }
283         debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9)
284
285         if(mpGiTree != NULL) delete mpGiTree;
286         char treeFlag = (1<<(this->mLbmInitId+4));
287         mpGiTree = new ntlTree( 
288                         15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
289                         scene, treeFlag );
290 }
291
292 /*****************************************************************************/
293 /*! destroy tree etc. when geometry init done */
294 void LbmSolverInterface::freeGeoTree() {
295         if(mpGiTree != NULL) {
296                 delete mpGiTree;
297                 mpGiTree = NULL;
298         }
299 }
300
301
302 int globGeoInitDebug = 0;
303 int globGICPIProblems = 0;
304 /*****************************************************************************/
305 /*! check for a certain flag type at position org */
306 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance) {
307         // shift ve ctors to avoid rounding errors
308         org += ntlVec3Gfx(0.0001);
309         ntlVec3Gfx dir = ntlVec3Gfx(1.0, 0.0, 0.0);
310         OId = -1;
311         ntlRay ray(org, dir, 0, 1.0, mpGlob);
312         //int insCnt = 0;
313         bool done = false;
314         bool inside = false;
315         vector<int> giObjFirstHistSide;
316         giObjFirstHistSide.resize( mpGiObjects->size() );
317         for(size_t i=0; i<mGiObjInside.size(); i++) { 
318                 mGiObjInside[i] = 0; 
319                 mGiObjDistance[i] = -1.0; 
320                 mGiObjSecondDist[i] = -1.0; 
321                 giObjFirstHistSide[i] = 0; 
322         }
323         // if not inside, return distance to first hit
324         gfxReal firstHit=-1.0;
325         int     firstOId = -1;
326         if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org<<" f"<<flags<<" acc"<<mAccurateGeoinit);
327         
328         if(mAccurateGeoinit) {
329                 while(!done) {
330                         // find first inside intersection
331                         ntlTriangle *triIns = NULL;
332                         distance = -1.0;
333                         ntlVec3Gfx normal(0.0);
334                         mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
335                         if(triIns) {
336                                 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
337                                 LbmFloat orientation = dot(normal, dir);
338                                 OId = triIns->getObjectId();
339                                 if(orientation<=0.0) {
340                                         // outside hit
341                                         normal *= -1.0;
342                                         mGiObjInside[OId]++;
343                                         if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = 1;
344                                         if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
345                                 } else {
346                                         // inside hit
347                                         mGiObjInside[OId]++;
348                                         if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
349                                         if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation);
350                                         if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = -1;
351                                 }
352                                 norg += normal * getVecEpsilon();
353                                 ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
354                                 // remember first hit distance, in case we're not 
355                                 // inside anything
356                                 if(firstHit<0.0) {
357                                         firstHit = distance;
358                                         firstOId = OId;
359                                 }
360                         } else {
361                                 // no more intersections... return false
362                                 done = true;
363                         }
364                 }
365
366                 distance = -1.0;
367                 for(size_t i=0; i<mGiObjInside.size(); i++) {
368                         if(mGiObjInside[i]>0) {
369                                 bool mess = false;
370                                 if((mGiObjInside[i]%2)==1) {
371                                         if(giObjFirstHistSide[i] != -1) mess=true;
372                                 } else {
373                                         if(giObjFirstHistSide[i] !=  1) mess=true;
374                                 }
375                                 if(mess) {
376                                         //errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
377                                         globGICPIProblems++;
378                                         mGiObjInside[i]++; // believe first hit side...
379                                 }
380                         }
381                 }
382                 for(size_t i=0; i<mGiObjInside.size(); i++) {
383                         if(globGeoInitDebug) errMsg("CHIII","i"<<i<<" ins="<<mGiObjInside[i]<<" t"<<mGiObjDistance[i]<<" d"<<distance);
384                         if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
385                                 if(  (distance<0.0)                             || // first intersection -> good
386                                           ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
387                                                 ) {                                             
388                                         distance = mGiObjDistance[i];
389                                         OId = i;
390                                         inside = true;
391                                 } 
392                         }
393                 }
394                 if(!inside) {
395                         distance = firstHit;
396                         OId = firstOId;
397                 }
398                 if(globGeoInitDebug) errMsg("CHIII","i"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
399
400                 return inside;
401         } else {
402
403                 // find first inside intersection
404                 ntlTriangle *triIns = NULL;
405                 distance = -1.0;
406                 ntlVec3Gfx normal(0.0);
407                 mpGiTree->intersectX(ray,distance,normal, triIns, flags, true);
408                 if(triIns) {
409                         // check outside intersect
410                         LbmFloat orientation = dot(normal, dir);
411                         if(orientation<=0.0) return false;
412
413                         OId = triIns->getObjectId();
414                         return true;
415                 }
416                 return false;
417         }
418 }
419 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse) {
420         // shift ve ctors to avoid rounding errors
421         org += ntlVec3Gfx(0.0001); //?
422         OId = -1;
423         ntlRay ray(org, dir, 0, 1.0, mpGlob);
424         //int insCnt = 0;
425         bool done = false;
426         bool inside = false;
427         for(size_t i=0; i<mGiObjInside.size(); i++) { 
428                 mGiObjInside[i] = 0; 
429                 mGiObjDistance[i] = -1.0; 
430                 mGiObjSecondDist[i] = -1.0; 
431         }
432         // if not inside, return distance to first hit
433         gfxReal firstHit=-1.0;
434         int     firstOId = -1;
435         thinHit = false;
436         
437         if(mAccurateGeoinit) {
438                 while(!done) {
439                         // find first inside intersection
440                         ntlTriangle *triIns = NULL;
441                         distance = -1.0;
442                         ntlVec3Gfx normal(0.0);
443                         mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
444                         if(triIns) {
445                                 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
446                                 LbmFloat orientation = dot(normal, dir);
447                                 OId = triIns->getObjectId();
448                                 if(orientation<=0.0) {
449                                         // outside hit
450                                         normal *= -1.0;
451                                         //mGiObjDistance[OId] = -1.0;
452                                         //errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
453                                 } else {
454                                         // inside hit
455                                         //if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
456                                         //errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
457                                         if(mGiObjInside[OId]==1) {
458                                                 // second inside hit
459                                                 if(mGiObjSecondDist[OId]<0.0) mGiObjSecondDist[OId] = distance;
460                                         }
461                                 }
462                                 mGiObjInside[OId]++;
463                                 // always store first hit for thin obj init
464                                 if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
465
466                                 norg += normal * getVecEpsilon();
467                                 ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
468                                 // remember first hit distance, in case we're not 
469                                 // inside anything
470                                 if(firstHit<0.0) {
471                                         firstHit = distance;
472                                         firstOId = OId;
473                                 }
474                         } else {
475                                 // no more intersections... return false
476                                 done = true;
477                                 //if(insCnt%2) inside=true;
478                         }
479                 }
480
481                 distance = -1.0;
482                 // standard inside check
483                 for(size_t i=0; i<mGiObjInside.size(); i++) {
484                         if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
485                                 if(  (distance<0.0)                             || // first intersection -> good
486                                           ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
487                                                 ) {                                             
488                                         distance = mGiObjDistance[i];
489                                         OId = i;
490                                         inside = true;
491                                 } 
492                         }
493                 }
494                 // now check for thin hits
495                 if(!inside) {
496                         distance = -1.0;
497                         for(size_t i=0; i<mGiObjInside.size(); i++) {
498                                 if((mGiObjInside[i]>=2)&&(mGiObjDistance[i]>0.0)&&(mGiObjSecondDist[i]>0.0)&&
499                                          (mGiObjDistance[i]<1.0*halfCellsize)&&(mGiObjSecondDist[i]<2.0*halfCellsize) ) {
500                                         if(  (distance<0.0)                             || // first intersection -> good
501                                                         ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
502                                                         ) {                                             
503                                                 distance = mGiObjDistance[i];
504                                                 OId = i;
505                                                 inside = true;
506                                                 thinHit = true;
507                                         } 
508                                 }
509                         }
510                 }
511                 if(!inside) {
512                         // check for hit in this cell, opposite to current dir (only recurse once)
513                         if(recurse) {
514                                 gfxReal r_distance;
515                                 int r_OId;
516                                 bool ret = geoInitCheckPointInside(org, dir*-1.0, flags, r_OId, r_distance, halfCellsize, thinHit, false);
517                                 if((ret)&&(thinHit)) {
518                                         OId = r_OId;
519                                         distance = 0.0; 
520                                         return true;
521                                 }
522                         }
523                 }
524                 // really no hit...
525                 if(!inside) {
526                         distance = firstHit;
527                         OId = firstOId;
528                         /*if((mGiObjDistance[OId]>0.0)&&(mGiObjSecondDist[OId]>0.0)) {
529                                 const gfxReal thisdist = mGiObjSecondDist[OId]-mGiObjDistance[OId];
530                                 // dont walk over this cell...
531                                 if(thisdist<halfCellsize) distance-=2.0*halfCellsize;
532                         } // ? */
533                 }
534                 //errMsg("CHIII","i"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
535
536                 return inside;
537         } else {
538
539                 // find first inside intersection
540                 ntlTriangle *triIns = NULL;
541                 distance = -1.0;
542                 ntlVec3Gfx normal(0.0);
543                 mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
544                 if(triIns) {
545                         // check outside intersect
546                         LbmFloat orientation = dot(normal, dir);
547                         if(orientation<=0.0) return false;
548
549                         OId = triIns->getObjectId();
550                         return true;
551                 }
552                 return false;
553         }
554 }
555
556
557 /*****************************************************************************/
558 ntlVec3Gfx LbmSolverInterface::getGeoMaxMovementVelocity(LbmFloat simtime, LbmFloat stepsize) {
559         ntlVec3Gfx max(0.0);
560         if(mpGlob == NULL) return max;
561         // mpGiObjects has to be inited here...
562         
563         for(int i=0; i< (int)mpGiObjects->size(); i++) {
564                 //errMsg("LbmSolverInterface::getGeoMaxMovementVelocity","i="<<i<<" "<< (*mpGiObjects)[i]->getName() ); // DEBUG
565                 if( (*mpGiObjects)[i]->getGeoInitType() & (FGI_FLUID|FGI_MBNDINFLOW) ){
566                         //ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
567                         ntlVec3Gfx orgvel = (*mpGiObjects)[i]->calculateMaxVel( simtime, simtime+stepsize );
568                         if( normNoSqrt(orgvel) > normNoSqrt(max) ) { max = orgvel; } 
569                         //errMsg("MVT","i"<<i<<" a"<< (*mpGiObjects)[i]->getInitialVelocity(simtime)<<" o"<<orgvel ); // DEBUG
570                         // TODO check if inflow and simtime 
571                         ntlVec3Gfx inivel = (*mpGiObjects)[i]->getInitialVelocity(simtime);
572                         if( normNoSqrt(inivel) > normNoSqrt(max) ) { max = inivel; } 
573                 }
574         }
575         errMsg("LbmSolverInterface::getGeoMaxMovementVelocity", "max="<< max ); // DEBUG
576         return max;
577 }
578
579
580
581
582 /*******************************************************************************/
583 /*! cell iteration functions */
584 /*******************************************************************************/
585
586
587
588
589 /*****************************************************************************/
590 //! add cell to mMarkedCells list
591 void 
592 LbmSolverInterface::addCellToMarkedList( CellIdentifierInterface *cid ) {
593         for(size_t i=0; i<mMarkedCells.size(); i++) {
594                 // check if cids alreay in
595                 if( mMarkedCells[i]->equal(cid) ) return;
596                 //mMarkedCells[i]->setEnd(false);
597         }
598         mMarkedCells.push_back( cid );
599         //cid->setEnd(true);
600 }
601
602 /*****************************************************************************/
603 //! marked cell iteration methods
604 CellIdentifierInterface* 
605 LbmSolverInterface::markedGetFirstCell( ) {
606         if(mMarkedCells.size() > 0){ return mMarkedCells[0]; }
607         return NULL;
608 }
609
610 CellIdentifierInterface* 
611 LbmSolverInterface::markedAdvanceCell() {
612         mMarkedCellIndex++;
613         if(mMarkedCellIndex>=(int)mMarkedCells.size()) return NULL;
614         return mMarkedCells[mMarkedCellIndex];
615 }
616
617 void LbmSolverInterface::markedClearList() {
618         // FIXME free cids?
619         mMarkedCells.clear();
620 }
621
622 /*******************************************************************************/
623 /*! string helper functions */
624 /*******************************************************************************/
625
626
627
628 // 32k
629 string convertSingleFlag2String(CellFlagType cflag) {
630         CellFlagType flag = cflag;
631         if(flag == CFUnused         ) return string("cCFUnused");
632         if(flag == CFEmpty          ) return string("cCFEmpty");      
633         if(flag == CFBnd            ) return string("cCFBnd");        
634         if(flag == CFBndNoslip      ) return string("cCFBndNoSlip");        
635         if(flag == CFBndFreeslip    ) return string("cCFBndFreeSlip");        
636         if(flag == CFBndPartslip    ) return string("cCFBndPartSlip");        
637         if(flag == CFNoInterpolSrc  ) return string("cCFNoInterpolSrc");
638         if(flag == CFFluid          ) return string("cCFFluid");      
639         if(flag == CFInter          ) return string("cCFInter");      
640         if(flag == CFNoNbFluid      ) return string("cCFNoNbFluid");  
641         if(flag == CFNoNbEmpty      ) return string("cCFNoNbEmpty");  
642         if(flag == CFNoDelete       ) return string("cCFNoDelete");   
643         if(flag == CFNoBndFluid     ) return string("cCFNoBndFluid"); 
644         if(flag == CFGrNorm         ) return string("cCFGrNorm");     
645         if(flag == CFGrFromFine     ) return string("cCFGrFromFine"); 
646         if(flag == CFGrFromCoarse   ) return string("cCFGrFromCoarse");
647         if(flag == CFGrCoarseInited ) return string("cCFGrCoarseInited");
648         if(flag == CFMbndInflow )     return string("cCFMbndInflow");
649         if(flag == CFMbndOutflow )    return string("cCFMbndOutflow");
650         if(flag == CFInvalid        ) return string("cfINVALID");
651
652         std::ostringstream mult;
653         int val = 0;
654         if(flag != 0) {
655                 while(! (flag&1) ) {
656                         flag = flag>>1;
657                         val++;
658                 }
659         } else {
660                 val = -1;
661         }
662         if(val>=24) {
663                 mult << "cfOID_" << (flag>>24) <<"_TYPE";
664         } else {
665                 mult << "cfUNKNOWN_" << val <<"_TYPE";
666         }
667         return mult.str();
668 }
669         
670 //! helper function to convert flag to string (for debuggin)
671 string convertCellFlagType2String( CellFlagType cflag ) {
672         int flag = cflag;
673
674         const int jmax = sizeof(CellFlagType)*8;
675         bool somefound = false;
676         std::ostringstream mult;
677         mult << "[";
678         for(int j=0; j<jmax ; j++) {
679                 if(flag& (1<<j)) {
680                         if(somefound) mult << "|";
681                         mult << j<<"<"<< convertSingleFlag2String( (CellFlagType)(1<<j) ); // this call should always be _non_-recursive
682                         somefound = true;
683                 }
684         };
685         mult << "]";
686
687         // return concatenated string
688         if(somefound) return mult.str();
689
690         // empty?
691         return string("[emptyCFT]");
692 }
693