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