initial commit of the fluid simulator.
[blender.git] / intern / elbeem / intern / lbminterface.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 "lbmdimensions.h" 
15 #include "lbminterface.h" 
16 #include "lbmfunctions.h" 
17 #include "ntl_scene.h"
18 #include "ntl_ray.h"
19 #include "typeslbm.h"
20
21
22 /*****************************************************************************/
23 //! common variables 
24
25 /*****************************************************************************/
26 /*! class for solver templating - 3D implementation D3Q19 */
27
28         //! how many dimensions?
29         const int LbmD3Q19::cDimension = 3;
30
31         // Wi factors for collide step 
32         const LbmFloat LbmD3Q19::cCollenZero    = (1.0/3.0);
33         const LbmFloat LbmD3Q19::cCollenOne     = (1.0/18.0);
34         const LbmFloat LbmD3Q19::cCollenSqrtTwo = (1.0/36.0);
35
36         //! threshold value for filled/emptied cells 
37         const LbmFloat LbmD3Q19::cMagicNr2    = 1.0005;
38         const LbmFloat LbmD3Q19::cMagicNr2Neg = -0.0005;
39         const LbmFloat LbmD3Q19::cMagicNr     = 1.010001;
40         const LbmFloat LbmD3Q19::cMagicNrNeg  = -0.010001;
41
42         //! size of a single set of distribution functions 
43         const int    LbmD3Q19::cDfNum      = 19;
44         //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
45         const int    LbmD3Q19::cDirNum     = 27;
46
47         //const string LbmD3Q19::dfString[ cDfNum ] = { 
48         const char* LbmD3Q19::dfString[ cDfNum ] = { 
49                 " C", " N"," S"," E"," W"," T"," B",
50                 "NE","NW","SE","SW",
51                 "NT","NB","ST","SB",
52                 "ET","EB","WT","WB"
53         };
54
55         const LbmD3Q19::dfDir LbmD3Q19::dfNorm[ cDfNum ] = { 
56                 cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB, 
57                 cDirNE, cDirNW, cDirSE, cDirSW, 
58                 cDirNT, cDirNB, cDirST, cDirSB, 
59                 cDirET, cDirEB, cDirWT, cDirWB
60         };
61
62         const LbmD3Q19::dfDir LbmD3Q19::dfInv[ cDfNum ] = { 
63                 cDirC,  cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
64                 cDirSW, cDirSE, cDirNW, cDirNE,
65                 cDirSB, cDirST, cDirNB, cDirNT, 
66                 cDirWB, cDirWT, cDirEB, cDirET
67         };
68
69         const int LbmD3Q19::dfRefX[ cDfNum ] = { 
70                 0,  0, 0, 0, 0, 0, 0,
71                 cDirSE, cDirSW, cDirNE, cDirNW,
72                 0, 0, 0, 0, 
73                 cDirEB, cDirET, cDirWB, cDirWT
74         };
75
76         const int LbmD3Q19::dfRefY[ cDfNum ] = { 
77                 0,  0, 0, 0, 0, 0, 0,
78                 cDirNW, cDirNE, cDirSW, cDirSE,
79                 cDirNB, cDirNT, cDirSB, cDirST,
80                 0, 0, 0, 0
81         };
82
83         const int LbmD3Q19::dfRefZ[ cDfNum ] = { 
84                 0,  0, 0, 0, 0, 0, 0,
85                 0, 0, 0, 0, 
86                 cDirST, cDirSB, cDirNT, cDirNB,
87                 cDirWT, cDirWB, cDirET, cDirEB
88         };
89
90         // Vector Order 3D:
91         //  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
92         //  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
93         //  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
94         //  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1
95
96         const int LbmD3Q19::dfVecX[ cDirNum ] = { 
97                 0, 0,0, 1,-1, 0,0,
98                 1,-1,1,-1,
99                 0,0,0,0,
100                 1,1,-1,-1,
101                  1,-1, 1,-1,
102                  1,-1, 1,-1,
103         };
104         const int LbmD3Q19::dfVecY[ cDirNum ] = { 
105                 0, 1,-1, 0,0,0,0,
106                 1,1,-1,-1,
107                 1,1,-1,-1,
108                 0,0,0,0,
109                  1, 1,-1,-1,
110                  1, 1,-1,-1
111         };
112         const int LbmD3Q19::dfVecZ[ cDirNum ] = { 
113                 0, 0,0,0,0,1,-1,
114                 0,0,0,0,
115                 1,-1,1,-1,
116                 1,-1,1,-1,
117                  1, 1, 1, 1,
118                 -1,-1,-1,-1
119         };
120
121         const LbmFloat LbmD3Q19::dfDvecX[ cDirNum ] = {
122                 0, 0,0, 1,-1, 0,0,
123                 1,-1,1,-1,
124                 0,0,0,0,
125                 1,1,-1,-1,
126                  1,-1, 1,-1,
127                  1,-1, 1,-1
128         };
129         const LbmFloat LbmD3Q19::dfDvecY[ cDirNum ] = {
130                 0, 1,-1, 0,0,0,0,
131                 1,1,-1,-1,
132                 1,1,-1,-1,
133                 0,0,0,0,
134                  1, 1,-1,-1,
135                  1, 1,-1,-1
136         };
137         const LbmFloat LbmD3Q19::dfDvecZ[ cDirNum ] = {
138                 0, 0,0,0,0,1,-1,
139                 0,0,0,0,
140                 1,-1,1,-1,
141                 1,-1,1,-1,
142                  1, 1, 1, 1,
143                 -1,-1,-1,-1
144         };
145
146         /* principal directions */
147         const int LbmD3Q19::princDirX[ 2*LbmD3Q19::cDimension ] = { 
148                 1,-1, 0,0, 0,0
149         };
150         const int LbmD3Q19::princDirY[ 2*LbmD3Q19::cDimension ] = { 
151                 0,0, 1,-1, 0,0
152         };
153         const int LbmD3Q19::princDirZ[ 2*LbmD3Q19::cDimension ] = { 
154                 0,0, 0,0, 1,-1
155         };
156
157         /*! arrays for les model coefficients, inited in lbmsolver constructor */
158         LbmFloat LbmD3Q19::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
159         LbmFloat LbmD3Q19::lesCoeffOffdiag[ cDimension ][ cDirNum ];
160
161
162         const LbmFloat LbmD3Q19::dfLength[ cDfNum ]= { 
163                 cCollenZero,
164                 cCollenOne, cCollenOne, cCollenOne, 
165                 cCollenOne, cCollenOne, cCollenOne,
166                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
167                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
168                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
169         };
170
171         /* precalculated equilibrium dfs, inited in lbmsolver constructor */
172         LbmFloat LbmD3Q19::dfEquil[ cDfNum ];
173
174 // D3Q19 end
175
176
177
178 /*****************************************************************************/
179 /*! class for solver templating - 2D implementation D2Q9 */
180
181                 //! how many dimensions?
182                 const int LbmD2Q9::cDimension = 2;
183
184                 //! Wi factors for collide step 
185                 const LbmFloat LbmD2Q9::cCollenZero    = (4.0/9.0);
186                 const LbmFloat LbmD2Q9::cCollenOne     = (1.0/9.0);
187                 const LbmFloat LbmD2Q9::cCollenSqrtTwo = (1.0/36.0);
188
189                 //! threshold value for filled/emptied cells 
190                 const LbmFloat LbmD2Q9::cMagicNr2    = 1.0005;
191                 const LbmFloat LbmD2Q9::cMagicNr2Neg = -0.0005;
192                 const LbmFloat LbmD2Q9::cMagicNr     = 1.010001;
193                 const LbmFloat LbmD2Q9::cMagicNrNeg  = -0.010001;
194
195                 //! size of a single set of distribution functions 
196                 const int LbmD2Q9::cDfNum  = 9;
197                 const int LbmD2Q9::cDirNum = 9;
198
199         //const string LbmD2Q9::dfString[ cDfNum ] = { 
200         const char* LbmD2Q9::dfString[ cDfNum ] = { 
201                 " C", 
202                 " N",   " S", " E", " W",
203                 "NE", "NW", "SE","SW" 
204         };
205
206         const LbmD2Q9::dfDir LbmD2Q9::dfNorm[ cDfNum ] = { 
207                 cDirC, 
208                 cDirN,  cDirS,  cDirE,  cDirW,
209                 cDirNE, cDirNW, cDirSE, cDirSW 
210         };
211
212         const LbmD2Q9::dfDir LbmD2Q9::dfInv[ cDfNum ] = { 
213                 cDirC,  
214                 cDirS,  cDirN,  cDirW,  cDirE,
215                 cDirSW, cDirSE, cDirNW, cDirNE 
216         };
217
218         const int LbmD2Q9::dfRefX[ cDfNum ] = { 
219                 0,  
220                 0,  0,  0,  0,
221                 cDirSE, cDirSW, cDirNE, cDirNW 
222         };
223
224         const int LbmD2Q9::dfRefY[ cDfNum ] = { 
225                 0,  
226                 0,  0,  0,  0,
227                 cDirNW, cDirNE, cDirSW, cDirSE 
228         };
229
230         const int LbmD2Q9::dfRefZ[ cDfNum ] = { 
231                 0,  0, 0, 0, 0,
232                 0, 0, 0, 0
233         };
234
235         // Vector Order 2D:
236         // 0  1 2  3  4  5  6 7  8
237         // 0, 0,0, 1,-1, 1,-1,1,-1 
238         // 0, 1,-1, 0,0, 1,1,-1,-1 
239
240         const int LbmD2Q9::dfVecX[ cDirNum ] = { 
241                 0, 
242                 0,0, 1,-1,
243                 1,-1,1,-1 
244         };
245         const int LbmD2Q9::dfVecY[ cDirNum ] = { 
246                 0, 
247                 1,-1, 0,0,
248                 1,1,-1,-1 
249         };
250         const int LbmD2Q9::dfVecZ[ cDirNum ] = { 
251                 0, 0,0,0,0, 0,0,0,0 
252         };
253
254         const LbmFloat LbmD2Q9::dfDvecX[ cDirNum ] = {
255                 0, 
256                 0,0, 1,-1,
257                 1,-1,1,-1 
258         };
259         const LbmFloat LbmD2Q9::dfDvecY[ cDirNum ] = {
260                 0, 
261                 1,-1, 0,0,
262                 1,1,-1,-1 
263         };
264         const LbmFloat LbmD2Q9::dfDvecZ[ cDirNum ] = {
265                 0, 0,0,0,0, 0,0,0,0 
266         };
267
268         const int LbmD2Q9::princDirX[ 2*LbmD2Q9::cDimension ] = { 
269                 1,-1, 0,0
270         };
271         const int LbmD2Q9::princDirY[ 2*LbmD2Q9::cDimension ] = { 
272                 0,0, 1,-1
273         };
274         const int LbmD2Q9::princDirZ[ 2*LbmD2Q9::cDimension ] = { 
275                 0,0, 0,0
276         };
277
278
279         /*! arrays for les model coefficients, inited in lbmsolver constructor */
280         LbmFloat LbmD2Q9::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
281         LbmFloat LbmD2Q9::lesCoeffOffdiag[ cDimension ][ cDirNum ];
282
283
284         const LbmFloat LbmD2Q9::dfLength[ cDfNum ]= { 
285                 cCollenZero,
286                 cCollenOne, cCollenOne, cCollenOne, cCollenOne, 
287                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
288         };
289
290         /* precalculated equilibrium dfs, inited in lbmsolver constructor */
291         LbmFloat LbmD2Q9::dfEquil[ cDfNum ];
292
293 // D2Q9 end
294
295
296
297 /******************************************************************************
298  * Interface Constructor
299  *****************************************************************************/
300 LbmSolverInterface::LbmSolverInterface() :
301         mPanic( false ),
302   mSizex(10), mSizey(10), mSizez(10), 
303   mStepCnt( 0 ),
304         mFixMass( 0.0 ),
305   mOmega( 1.0 ),
306   mGravity(0.0),
307         mSurfaceTension( 0.0 ), 
308   mInitialMass (0.0), 
309         mBoundaryEast(  (CellFlagType)(CFBnd) ),mBoundaryWest( (CellFlagType)(CFBnd) ),mBoundaryNorth(  (CellFlagType)(CFBnd) ),
310         mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop(  (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
311   mInitDone( false ),
312   mInitDensityGradient( false ),
313         mpAttrs( NULL ), mpParam( NULL ),
314         mNumParticlesLost(0), mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
315         mDebugVelScale( 1.0 ), mNodeInfoString("+"),
316         mRandom( 5123 ),
317         mvGeoStart(-1.0), mvGeoEnd(1.0),
318         mPerformGeoInit( false ),
319         mName("lbm_default") ,
320         mpIso( NULL ), mIsoValue(0.49999),
321         mSilent(false) , 
322         mGeoInitId( 0 ),
323         mpGiTree( NULL ),
324         mAccurateGeoinit(0),
325         mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL )
326 {
327 #if ELBEEM_BLENDER==1
328         mSilent = true;
329 #endif
330 }
331
332
333 /*******************************************************************************/
334 /*! parse a boundary flag string */
335 CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) {
336         string val  = mpAttrs->readString(name, "", source, target, needed);
337         if(!strcmp(val.c_str(),"")) {
338                 // no value given...
339                 return CFEmpty;
340         }
341         if(!strcmp(val.c_str(),"bnd_no")) {
342                 return (CellFlagType)( CFBnd );
343         }
344         if(!strcmp(val.c_str(),"bnd_free")) {
345                 return (CellFlagType)( CFBnd );
346         }
347         if(!strcmp(val.c_str(),"fluid")) {
348                 /* might be used for some in/out flow cases */
349                 return (CellFlagType)( CFFluid );
350         }
351         errorOut("LbmStdSolver::readBoundaryFlagInt error: Invalid value '"<<val<<"' " );
352 # if LBM_STRICT_DEBUG==1
353         exit(1);
354 # endif
355         return defaultValue;
356 }
357
358 /*******************************************************************************/
359 /*! parse standard attributes */
360 void LbmSolverInterface::parseStdAttrList() {
361         if(!mpAttrs) {
362                 errorOut("LbmStdSolver::parseAttrList error: mpAttrs pointer not initialized!");
363                 exit(1); }
364
365         // st currently unused
366         //mSurfaceTension  = mpAttrs->readFloat("d_surfacetension",  mSurfaceTension, "LbmStdSolver", "mSurfaceTension", false);
367         mBoundaryEast  = readBoundaryFlagInt("boundary_east",  mBoundaryEast, "LbmStdSolver", "mBoundaryEast", false);
368         mBoundaryWest  = readBoundaryFlagInt("boundary_west",  mBoundaryWest, "LbmStdSolver", "mBoundaryWest", false);
369         mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmStdSolver", "mBoundaryNorth", false);
370         mBoundarySouth = readBoundaryFlagInt("boundary_south", mBoundarySouth,"LbmStdSolver", "mBoundarySouth", false);
371         mBoundaryTop   = readBoundaryFlagInt("boundary_top",   mBoundaryTop,"LbmStdSolver", "mBoundaryTop", false);
372         mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmStdSolver", "mBoundaryBottom", false);
373
374         LbmVec sizeVec(mSizex,mSizey,mSizez);
375         sizeVec = vec2L( mpAttrs->readVec3d("size",  vec2P(sizeVec), "LbmStdSolver", "sizeVec", false) );
376         mSizex = (int)sizeVec[0]; 
377         mSizey = (int)sizeVec[1]; 
378         mSizez = (int)sizeVec[2];
379         mpParam->setSize(mSizex, mSizey, mSizez ); // param needs size in any case
380
381         mInitDensityGradient = mpAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmStdSolver", "mInitDensityGradient", false);
382         mPerformGeoInit = mpAttrs->readBool("geoinit", mPerformGeoInit,"LbmStdSolver", "mPerformGeoInit", false);
383         mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"LbmStdSolver", "mGeoInitId", false);
384         mIsoValue = mpAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false );
385
386         mDebugVelScale = mpAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmStdSolver", "mDebugVelScale", false);
387         mNodeInfoString = mpAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false );
388 }
389
390
391 /*******************************************************************************/
392 /*! geometry initialization */
393 /*******************************************************************************/
394
395 /*****************************************************************************/
396 /*! init tree for certain geometry init */
397 void LbmSolverInterface::initGeoTree(int id) {
398         if(mpGlob == NULL) { errorOut("LbmSolverInterface::initGeoTree error: Requires globals!"); exit(1); }
399         mGeoInitId = id;
400         ntlScene *scene = mpGlob->getScene();
401         mpGiObjects = scene->getObjects();
402         mGiObjInside.resize( mpGiObjects->size() );
403         mGiObjDistance.resize( mpGiObjects->size() );
404         for(size_t i=0; i<mpGiObjects->size(); i++) { 
405                 if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
406         }
407         debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9)
408
409         if(mpGiTree != NULL) delete mpGiTree;
410         char treeFlag = (1<<(mGeoInitId+4));
411         mpGiTree = new ntlTree( 20, 4, // warning - fixed values for depth & maxtriangles here...
412                                                                                                 scene, treeFlag );
413 }
414
415 /*****************************************************************************/
416 /*! destroy tree etc. when geometry init done */
417 void LbmSolverInterface::freeGeoTree() {
418         if(mpGiTree != NULL) {
419                 delete mpGiTree;
420                 mpGiTree = NULL;
421         }
422 }
423
424
425 /*****************************************************************************/
426 /*! check for a certain flag type at position org */
427 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance) {
428         // shift ve ctors to avoid rounding errors
429         org += ntlVec3Gfx(0.0001);
430         ntlVec3Gfx dir = ntlVec3Gfx(0.999999,0.0,0.0);
431         OId = -1;
432         ntlRay ray(org, dir, 0, 1.0, mpGlob);
433         //int insCnt = 0;
434         bool done = false;
435         bool inside = false;
436         //errMsg("III"," start org"<<org<<" dir"<<dir);
437         //int insID = ray.getID();
438         for(size_t i=0; i<mGiObjInside.size(); i++) { mGiObjInside[i] = 0; }
439         // if not inside, return distance to first hit
440         gfxReal firstHit=-1.0;
441         
442         if(mAccurateGeoinit) {
443                 while(!done) {
444                         // find first inside intersection
445                         ntlTriangle *triIns = NULL;
446                         distance = -1.0;
447                         ntlVec3Gfx normal(0.0);
448                         mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
449                         if(triIns) {
450                                 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
451                                 LbmFloat orientation = dot(normal, dir);
452                                 OId = triIns->getObjectId();
453                                 if(orientation<=0.0) {
454                                         // outside hit
455                                         normal *= -1.0;
456                                         mGiObjInside[OId]++;
457                                         mGiObjDistance[OId] = -1.0;
458                                         //errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
459                                 } else {
460                                         // inside hit
461                                         mGiObjInside[OId]++;
462                                         mGiObjDistance[OId] = distance;
463                                         //errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
464                                 }
465                                 norg += normal * getVecEpsilon();
466                                 ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
467                                 if(firstHit<0.0) firstHit = distance;
468                                 //if((OId<0) && ())
469                                 //insCnt++;
470                                 /*
471                                 // check outside intersect
472                                 LbmFloat orientation = dot(normal, dir);
473                                 if(orientation<=0.0) {
474                                 // do more thorough checks... advance ray
475                                 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
476                                 norg += normal * (-1.0 * getVecEpsilon());
477                                 ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
478                                 insCnt++;
479                                 errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" insCnt"<<insCnt);
480                                 } else {
481                                 if(insCnt>0) {
482                                 // we have entered this obj before?
483                                 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
484                                 norg += normal * (-1.0 * getVecEpsilon());
485                                 ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
486                                 insCnt--;
487                                 errMsg("IIIS"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" insCnt"<<insCnt);
488                                 } else {
489                                 // first inside intersection -> ok
490                                 OId = triIns->getObjectId();
491                                 done = inside = true;
492                                 }
493                                 }
494                                 */
495                         } else {
496                                 // no more intersections... return false
497                                 done = true;
498                                 //if(insCnt%2) inside=true;
499                         }
500                 }
501
502                 distance = -1.0;
503                 for(size_t i=0; i<mGiObjInside.size(); i++) {
504                         //errMsg("CHIII","i"<<i<<" ins="<<mGiObjInside[i]<<" t"<<mGiObjDistance[i]<<" d"<<distance);
505                         if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
506                                 if(distance<0.0) {
507                                         // first intersection -> good
508                                         distance = mGiObjDistance[i];
509                                         OId = i;
510                                         inside = true;
511                                 } else {
512                                         if(distance>mGiObjDistance[i]) {
513                                                 // more than one intersection -> use closest one
514                                                 distance = mGiObjDistance[i];
515                                                 OId = i;
516                                                 inside = true;
517                                         }
518                                 }
519                         }
520                 }
521                 if(!inside) {
522                         distance = firstHit;
523                 }
524
525                 return inside;
526         } else {
527
528                 // find first inside intersection
529                 ntlTriangle *triIns = NULL;
530                 distance = -1.0;
531                 ntlVec3Gfx normal(0.0);
532                 mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
533                 if(triIns) {
534                         // check outside intersect
535                         LbmFloat orientation = dot(normal, dir);
536                         if(orientation<=0.0) return false;
537
538                         OId = triIns->getObjectId();
539                         return true;
540                 }
541                 return false;
542         }
543 }
544
545 /*****************************************************************************/
546 /*! get max. velocity of all objects to initialize as fluid regions */
547 ntlVec3Gfx LbmSolverInterface::getGeoMaxInitialVelocity() {
548         ntlScene *scene = mpGlob->getScene();
549         mpGiObjects = scene->getObjects();
550         ntlVec3Gfx max(0.0);
551         
552         for(int i=0; i< (int)mpGiObjects->size(); i++) {
553                 if( (*mpGiObjects)[i]->getGeoInitType() & FGI_FLUID ){
554                         ntlVec3Gfx ovel = (*mpGiObjects)[i]->getInitialVelocity();
555                         if( normNoSqrt(ovel) > normNoSqrt(max) ) { max = ovel; } 
556                         //errMsg("IVT","i"<<i<<" "<< (*mpGiObjects)[i]->getInitialVelocity() ); // DEBUG
557                 }
558         }
559         //errMsg("IVT","max "<<" "<< max ); // DEBUG
560         // unused !? mGiInsideCnt.resize( mpGiObjects->size() );
561
562         return max;
563 }
564
565
566 /*******************************************************************************/
567 /*! "traditional" initialization */
568 /*******************************************************************************/
569
570
571 /*****************************************************************************/
572 // handle generic test cases (currently only reset geo init)
573 bool LbmSolverInterface::initGenericTestCases() {
574         bool initTypeFound = false;
575         LbmSolverInterface::CellIdentifier cid = getFirstCell();
576         // deprecated! - only check for invalid cells...
577
578         // this is the default init - check if the geometry flag init didnt
579         initTypeFound = true;
580
581         while(noEndCell(cid)) {
582                 // check node
583                 if( (getCellFlag(cid,0)==CFInvalid) || (getCellFlag(cid,1)==CFInvalid) ) {
584                         warnMsg("LbmSolverInterface::initGenericTestCases","GeoInit produced invalid Flag at "<<cid->getAsString()<<"!" );
585                 }
586                 advanceCell( cid );
587         }
588
589         deleteCellIterator( &cid );
590         return initTypeFound;
591 }
592
593
594 /*******************************************************************************/
595 /*! cell iteration functions */
596 /*******************************************************************************/
597
598
599
600
601 /*****************************************************************************/
602 //! add cell to mMarkedCells list
603 void 
604 LbmSolverInterface::addCellToMarkedList( CellIdentifierInterface *cid ) {
605         for(size_t i=0; i<mMarkedCells.size(); i++) {
606                 // check if cids alreay in
607                 if( mMarkedCells[i]->equal(cid) ) return;
608                 mMarkedCells[i]->setEnd(false);
609         }
610         mMarkedCells.push_back( cid );
611         cid->setEnd(true);
612 }
613
614 /*****************************************************************************/
615 //! marked cell iteration methods
616 CellIdentifierInterface* 
617 LbmSolverInterface::markedGetFirstCell( ) {
618         /*MarkedCellIdentifier *newcid = new MarkedCellIdentifier();
619         if(mMarkedCells.size() < 1){ 
620                 newcid->setEnd( true );
621         }       else {
622                 newcid->mpCell = mMarkedCells[0];
623         }
624         return newcid;*/
625         return NULL;
626 }
627
628 void 
629 LbmSolverInterface::markedAdvanceCell( CellIdentifierInterface* basecid ) {
630         if(!basecid) return;
631         basecid->setEnd( true );
632         /*MarkedCellIdentifier *cid = dynamic_cast<MarkedCellIdentifier*>( basecid );
633         CellIdentifierInterface *cid = basecid;
634         cid->mIndex++;
635         if(cid->mIndex >= (int)mMarkedCells.size()) {
636                 cid->setEnd( true );
637         }
638         cid->mpCell = mMarkedCells[ cid->mIndex ];
639         */
640 }
641
642 bool 
643 LbmSolverInterface::markedNoEndCell( CellIdentifierInterface* cid ) {
644         if(!cid) return false;
645         return(! cid->getEnd() );
646 }
647
648 void LbmSolverInterface::markedClearList() {
649         // FIXME free cids?
650         mMarkedCells.clear();
651 }
652
653 /*******************************************************************************/
654 /*! string helper functions */
655 /*******************************************************************************/
656
657
658
659 // 32k
660 std::string convertSingleFlag2String(CellFlagType cflag) {
661         CellFlagType flag = cflag;
662         if(flag == CFUnused         ) return string("cCFUnused");
663         if(flag == CFEmpty          ) return string("cCFEmpty");      
664         if(flag == CFBnd            ) return string("cCFBnd");        
665         if(flag == CFNoInterpolSrc  ) return string("cCFNoInterpolSrc");
666         if(flag == CFFluid          ) return string("cCFFluid");      
667         if(flag == CFInter          ) return string("cCFInter");      
668         if(flag == CFNoNbFluid      ) return string("cCFNoNbFluid");  
669         if(flag == CFNoNbEmpty      ) return string("cCFNoNbEmpty");  
670         if(flag == CFNoDelete       ) return string("cCFNoDelete");   
671         if(flag == CFNoBndFluid     ) return string("cCFNoBndFluid"); 
672         if(flag == CFBndMARK        ) return string("cCFBndMARK");     
673         if(flag == CFGrNorm         ) return string("cCFGrNorm");     
674         if(flag == CFGrFromFine     ) return string("cCFGrFromFine"); 
675         if(flag == CFGrFromCoarse   ) return string("cCFGrFromCoarse");
676         if(flag == CFGrCoarseInited ) return string("cCFGrCoarseInited");
677         if(flag == CFInvalid        ) return string("cfINVALID");
678
679         std::ostringstream mult;
680         int val = 0;
681         if(flag != 0) {
682                 while(! (flag&1) ) {
683                         flag = flag>>1;
684                         val++;
685                 }
686         } else {
687                 val = -1;
688         }
689         mult << "cfUNKNOWN_" << val <<"_TYPE";
690         return mult.str();
691 }
692         
693 //! helper function to convert flag to string (for debuggin)
694 std::string convertCellFlagType2String( CellFlagType cflag ) {
695         int flag = cflag;
696
697         const int jmax = sizeof(CellFlagType)*8;
698         bool somefound = false;
699         std::ostringstream mult;
700         mult << "[";
701         for(int j=0; j<jmax ; j++) {
702                 if(flag& (1<<j)) {
703                         if(somefound) mult << "|";
704                         mult << convertSingleFlag2String( (CellFlagType)(1<<j) ); // this call should always be _non_-recursive
705                         somefound = true;
706                 }
707         };
708         mult << "]";
709
710         // return concatenated string
711         if(somefound) return mult.str();
712
713         // empty?
714         return string("[emptyCFT]");
715 }
716