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