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