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