|
Blender
V2.59
|
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method 00007 * All code distributed as part of El'Beem is covered by the version 2 of the 00008 * GNU General Public License. See the file COPYING for details. 00009 * Copyright 2003-2006 Nils Thuerey 00010 * 00011 * Combined 2D/3D Lattice Boltzmann Interface Class 00012 * contains stuff to be statically compiled 00013 * 00014 *****************************************************************************/ 00015 00016 /* LBM Files */ 00017 #include "ntl_matrices.h" 00018 #include "solver_interface.h" 00019 #include "ntl_ray.h" 00020 #include "ntl_world.h" 00021 #include "elbeem.h" 00022 00023 00024 00025 00026 /****************************************************************************** 00027 * Interface Constructor 00028 *****************************************************************************/ 00029 LbmSolverInterface::LbmSolverInterface() : 00030 mPanic( false ), 00031 mSizex(10), mSizey(10), mSizez(10), 00032 mAllfluid(false), mStepCnt( 0 ), 00033 mFixMass( 0.0 ), 00034 mOmega( 1.0 ), 00035 mGravity(0.0), 00036 mSurfaceTension( 0.0 ), 00037 mBoundaryEast( (CellFlagType)(CFBnd) ),mBoundaryWest( (CellFlagType)(CFBnd) ),mBoundaryNorth( (CellFlagType)(CFBnd) ), 00038 mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop( (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ), 00039 mInitDone( false ), 00040 mInitDensityGradient( false ), 00041 mpSifAttrs( NULL ), mpSifSwsAttrs(NULL), mpParam( NULL ), mpParticles(NULL), 00042 mNumParticlesLost(0), 00043 mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0), 00044 mDebugVelScale( 0.01 ), mNodeInfoString("+"), 00045 mvGeoStart(-1.0), mvGeoEnd(1.0), mpSimTrafo(NULL), 00046 mAccurateGeoinit(0), 00047 mName("lbm_default") , 00048 mpIso( NULL ), mIsoValue(0.499), 00049 mSilent(false) , 00050 mLbmInitId(1) , 00051 mpGiTree( NULL ), 00052 mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ), 00053 mRefinementDesired(0), 00054 mOutputSurfacePreview(0), mPreviewFactor(0.25), 00055 mSmoothSurface(1.0), mSmoothNormals(1.0), 00056 mIsoSubdivs(1), mPartGenProb(0.), 00057 mDumpVelocities(false), 00058 mMarkedCells(), mMarkedCellIndex(0), 00059 mDomainBound("noslip"), mDomainPartSlipValue(0.1), 00060 mFarFieldSize(0.), 00061 mPartDropMassSub(0.1), // good default 00062 mPartUsePhysModel(false), 00063 mTForceStrength(0.0), 00064 mCppfStage(0), 00065 mDumpRawText(false), 00066 mDumpRawBinary(false), 00067 mDumpRawBinaryZip(true) 00068 { 00069 #if ELBEEM_PLUGIN==1 00070 if(gDebugLevel<=1) setSilent(true); 00071 #endif 00072 mpSimTrafo = new ntlMat4Gfx(0.0); 00073 mpSimTrafo->initId(); 00074 00075 if(getenv("ELBEEM_RAWDEBUGDUMP")) { 00076 debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2); 00077 mDumpRawText = true; 00078 } 00079 00080 if(getenv("ELBEEM_BINDEBUGDUMP")) { 00081 debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2); 00082 mDumpRawBinary = true; 00083 } 00084 } 00085 00086 LbmSolverInterface::~LbmSolverInterface() 00087 { 00088 if(mpSimTrafo) delete mpSimTrafo; 00089 } 00090 00091 00092 /****************************************************************************** 00093 * initialize correct grid sizes given a geometric bounding box 00094 * and desired grid resolutions, all params except maxrefine 00095 * will be modified 00096 *****************************************************************************/ 00097 void initGridSizes(int &sizex, int &sizey, int &sizez, 00098 ntlVec3Gfx &geoStart, ntlVec3Gfx &geoEnd, 00099 int mMaxRefine, bool parallel) 00100 { 00101 // fix size inits to force cubic cells and mult4 level dimensions 00102 const int debugGridsizeInit = 1; 00103 if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" "<<geoStart<<","<<geoEnd ,10); 00104 00105 int maxGridSize = sizex; // get max size 00106 if(sizey>maxGridSize) maxGridSize = sizey; 00107 if(sizez>maxGridSize) maxGridSize = sizez; 00108 LbmFloat maxGeoSize = (geoEnd[0]-geoStart[0]); // get max size 00109 if((geoEnd[1]-geoStart[1])>maxGeoSize) maxGeoSize = (geoEnd[1]-geoStart[1]); 00110 if((geoEnd[2]-geoStart[2])>maxGeoSize) maxGeoSize = (geoEnd[2]-geoStart[2]); 00111 // FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though 00112 LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize); 00113 if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Start:"<<geoStart<<" End:"<<geoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10); 00114 // force grid sizes according to geom. size, rounded 00115 sizex = (int) ((geoEnd[0]-geoStart[0]) / cellSize +0.5); 00116 sizey = (int) ((geoEnd[1]-geoStart[1]) / cellSize +0.5); 00117 sizez = (int) ((geoEnd[2]-geoStart[2]) / cellSize +0.5); 00118 // match refinement sizes, round downwards to multiple of 4 00119 int sizeMask = 0; 00120 int maskBits = mMaxRefine; 00121 if(parallel==1) maskBits+=2; 00122 for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); } 00123 00124 // at least size 4 on coarsest level 00125 int minSize = 4<<mMaxRefine; //(maskBits+2); 00126 if(sizex<minSize) sizex = minSize; 00127 if(sizey<minSize) sizey = minSize; 00128 if(sizez<minSize) sizez = minSize; 00129 00130 sizeMask = ~sizeMask; 00131 if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask)<<", maskBits:"<<maskBits<<",minSize:"<<minSize ,10); 00132 sizex &= sizeMask; 00133 sizey &= sizeMask; 00134 sizez &= sizeMask; 00135 // debug 00136 int nextinc = (~sizeMask)+1; 00137 debMsgStd("initGridSizes",DM_MSG,"MPTeffMLtest, curr size:"<<PRINT_VEC(sizex,sizey,sizez)<<", next size:"<<PRINT_VEC(sizex+nextinc,sizey+nextinc,sizez+nextinc) ,10); 00138 00139 // force geom size to match rounded/modified grid sizes 00140 geoEnd[0] = geoStart[0] + cellSize*(LbmFloat)sizex; 00141 geoEnd[1] = geoStart[1] + cellSize*(LbmFloat)sizey; 00142 geoEnd[2] = geoStart[2] + cellSize*(LbmFloat)sizez; 00143 } 00144 00145 void calculateMemreqEstimate( int resx,int resy,int resz, 00146 int refine, float farfield, 00147 double *reqret, double *reqretFine, string *reqstr) { 00148 // debug estimation? 00149 const bool debugMemEst = true; 00150 // COMPRESSGRIDS define is not available here, make sure it matches 00151 const bool useGridComp = true; 00152 // make sure we can handle bid numbers here... all double 00153 double memCnt = 0.0; 00154 double ddTotalNum = (double)dTotalNum; 00155 if(reqretFine) *reqretFine = -1.; 00156 00157 double currResx = (double)resx; 00158 double currResy = (double)resy; 00159 double currResz = (double)resz; 00160 double rcellSize = ((currResx*currResy*currResz) *ddTotalNum); 00161 memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0); 00162 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"res:"<<PRINT_VEC(currResx,currResy,currResz)<<" rcellSize:"<<rcellSize<<" mc:"<<memCnt, 10); 00163 if(!useGridComp) { 00164 memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0); 00165 if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0); 00166 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," no-comp, mc:"<<memCnt, 10); 00167 } else { 00168 double compressOffset = (double)(currResx*currResy*ddTotalNum*2.0); 00169 memCnt += (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0)); 00170 if(reqretFine) *reqretFine = (double)(sizeof(LbmFloat) * (rcellSize+compressOffset +4.0)); 00171 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG," w-comp, mc:"<<memCnt, 10); 00172 } 00173 for(int i=refine-1; i>=0; i--) { 00174 currResx /= 2.0; 00175 currResy /= 2.0; 00176 currResz /= 2.0; 00177 rcellSize = ((currResz*currResy*currResx) *ddTotalNum); 00178 memCnt += (double)(sizeof(CellFlagType) * (rcellSize/ddTotalNum +4.0) *2.0); 00179 memCnt += (double)(sizeof(LbmFloat) * (rcellSize +4.0) *2.0); 00180 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"refine "<<i<<", mc:"<<memCnt, 10); 00181 } 00182 00183 // isosurface memory, use orig res values 00184 if(farfield>0.) { 00185 memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) ); 00186 } else { 00187 // ignore 3 int slices... 00188 memCnt += (double)( ( sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) ); 00189 } 00190 if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"iso, mc:"<<memCnt, 10); 00191 00192 // cpdata init check missing... 00193 00194 double memd = memCnt; 00195 const char *sizeStr = ""; 00196 const double sfac = 1024.0; 00197 if(memd>sfac){ memd /= sfac; sizeStr="KB"; } 00198 if(memd>sfac){ memd /= sfac; sizeStr="MB"; } 00199 if(memd>sfac){ memd /= sfac; sizeStr="GB"; } 00200 if(memd>sfac){ memd /= sfac; sizeStr="TB"; } 00201 00202 // return values 00203 std::ostringstream ret; 00204 if(memCnt< 1024.0*1024.0) { 00205 // show full MBs 00206 ret << (ceil(memd)); 00207 } else { 00208 // two digits for anything larger than MB 00209 ret << (ceil(memd*100.0)/100.0); 00210 } 00211 ret << " "<< sizeStr; 00212 *reqret = memCnt; 00213 *reqstr = ret.str(); 00214 //debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4); 00215 } 00216 00217 void LbmSolverInterface::initDomainTrafo(float *mat) { 00218 mpSimTrafo->initArrayCheck(mat); 00219 } 00220 00221 /*******************************************************************************/ 00223 CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) { 00224 string val = mpSifAttrs->readString(name, "", source, target, needed); 00225 if(!strcmp(val.c_str(),"")) { 00226 // no value given... 00227 return CFEmpty; 00228 } 00229 if(!strcmp(val.c_str(),"bnd_no")) { 00230 return (CellFlagType)( CFBnd ); 00231 } 00232 if(!strcmp(val.c_str(),"bnd_free")) { 00233 return (CellFlagType)( CFBnd ); 00234 } 00235 if(!strcmp(val.c_str(),"fluid")) { 00236 /* might be used for some in/out flow cases */ 00237 return (CellFlagType)( CFFluid ); 00238 } 00239 errMsg("LbmSolverInterface::readBoundaryFlagInt","Invalid value '"<<val<<"' " ); 00240 # if FSGR_STRICT_DEBUG==1 00241 errFatal("readBoundaryFlagInt","Strict abort..."<<val, SIMWORLD_INITERROR); 00242 # endif 00243 return defaultValue; 00244 } 00245 00246 /*******************************************************************************/ 00248 void LbmSolverInterface::parseStdAttrList() { 00249 if(!mpSifAttrs) { 00250 errFatal("LbmSolverInterface::parseAttrList","mpSifAttrs pointer not initialized!",SIMWORLD_INITERROR); 00251 return; } 00252 00253 // st currently unused 00254 mBoundaryEast = readBoundaryFlagInt("boundary_east", mBoundaryEast, "LbmSolverInterface", "mBoundaryEast", false); 00255 mBoundaryWest = readBoundaryFlagInt("boundary_west", mBoundaryWest, "LbmSolverInterface", "mBoundaryWest", false); 00256 mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmSolverInterface", "mBoundaryNorth", false); 00257 mBoundarySouth = readBoundaryFlagInt("boundary_south", mBoundarySouth,"LbmSolverInterface", "mBoundarySouth", false); 00258 mBoundaryTop = readBoundaryFlagInt("boundary_top", mBoundaryTop,"LbmSolverInterface", "mBoundaryTop", false); 00259 mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmSolverInterface", "mBoundaryBottom", false); 00260 00261 mpSifAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo ); 00262 00263 LbmVec sizeVec(mSizex,mSizey,mSizez); 00264 sizeVec = vec2L( mpSifAttrs->readVec3d("size", vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) ); 00265 mSizex = (int)sizeVec[0]; 00266 mSizey = (int)sizeVec[1]; 00267 mSizez = (int)sizeVec[2]; 00268 // param needs size in any case 00269 // test solvers might not have mpPara, though 00270 if(mpParam) mpParam->setSize(mSizex, mSizey, mSizez ); 00271 00272 mInitDensityGradient = mpSifAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false); 00273 mIsoValue = mpSifAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false ); 00274 00275 mDebugVelScale = mpSifAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false); 00276 mNodeInfoString = mpSifAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false ); 00277 00278 mDumpVelocities = mpSifAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false ); 00279 if(getenv("ELBEEM_DUMPVELOCITIES")) { 00280 int get = atoi(getenv("ELBEEM_DUMPVELOCITIES")); 00281 if((get==0)||(get==1)) { 00282 mDumpVelocities = get; 00283 debMsgStd("LbmSolverInterface::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_DUMPVELOCITIES to set mDumpVelocities to "<<get<<","<<mDumpVelocities,8); 00284 } 00285 } 00286 00287 // new test vars 00288 // mTForceStrength set from test solver 00289 mFarFieldSize = mpSifAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false); 00290 // old compat 00291 float sizeScale = mpSifAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false); 00292 if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); } 00293 00294 mPartDropMassSub = mpSifAttrs->readFloat("part_dropmass", mPartDropMassSub,"LbmSolverInterface", "mPartDropMassSub", false); 00295 00296 mCppfStage = mpSifAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false); 00297 mPartGenProb = mpSifAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false); 00298 mPartUsePhysModel = mpSifAttrs->readBool("partusephysmodel", mPartUsePhysModel,"LbmFsgrSolver", "mPartUsePhysModel", false); 00299 mIsoSubdivs = mpSifAttrs->readInt("isosubdivs", mIsoSubdivs,"LbmFsgrSolver", "mIsoSubdivs", false); 00300 } 00301 00302 00303 /*******************************************************************************/ 00305 /*******************************************************************************/ 00306 00307 /*****************************************************************************/ 00309 void LbmSolverInterface::initGeoTree() { 00310 if(mpGlob == NULL) { errFatal("LbmSolverInterface::initGeoTree","Requires globals!",SIMWORLD_INITERROR); return; } 00311 ntlScene *scene = mpGlob->getSimScene(); 00312 mpGiObjects = scene->getObjects(); 00313 mGiObjInside.resize( mpGiObjects->size() ); 00314 mGiObjDistance.resize( mpGiObjects->size() ); 00315 mGiObjSecondDist.resize( mpGiObjects->size() ); 00316 for(size_t i=0; i<mpGiObjects->size(); i++) { 00317 if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true; 00318 //debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 ); 00319 } 00320 debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9) 00321 debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Objects: "<<mpGiObjects->size() ,10); 00322 00323 if(mpGiTree != NULL) delete mpGiTree; 00324 char treeFlag = (1<<(this->mLbmInitId+4)); 00325 mpGiTree = new ntlTree( 00326 # if FSGR_STRICT_DEBUG!=1 00327 15, 8, // TREEwarning - fixed values for depth & maxtriangles here... 00328 # else // FSGR_STRICT_DEBUG==1 00329 10, 20, // reduced/slower debugging values 00330 # endif // FSGR_STRICT_DEBUG==1 00331 scene, treeFlag ); 00332 } 00333 00334 /*****************************************************************************/ 00336 void LbmSolverInterface::freeGeoTree() { 00337 if(mpGiTree != NULL) { 00338 delete mpGiTree; 00339 mpGiTree = NULL; 00340 } 00341 } 00342 00343 00344 int globGeoInitDebug = 0; 00345 int globGICPIProblems = 0; 00346 /*****************************************************************************/ 00348 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance, int shootDir) { 00349 // shift ve ctors to avoid rounding errors 00350 org += ntlVec3Gfx(0.0001); 00351 OId = -1; 00352 00353 // select shooting direction 00354 ntlVec3Gfx dir = ntlVec3Gfx(1., 0., 0.); 00355 if(shootDir==1) dir = ntlVec3Gfx(0., 1., 0.); 00356 else if(shootDir==2) dir = ntlVec3Gfx(0., 0., 1.); 00357 else if(shootDir==-2) dir = ntlVec3Gfx(0., 0., -1.); 00358 else if(shootDir==-1) dir = ntlVec3Gfx(0., -1., 0.); 00359 00360 ntlRay ray(org, dir, 0, 1.0, mpGlob); 00361 bool done = false; 00362 bool inside = false; 00363 00364 vector<int> giObjFirstHistSide; 00365 giObjFirstHistSide.resize( mpGiObjects->size() ); 00366 for(size_t i=0; i<mGiObjInside.size(); i++) { 00367 mGiObjInside[i] = 0; 00368 mGiObjDistance[i] = -1.0; 00369 mGiObjSecondDist[i] = -1.0; 00370 giObjFirstHistSide[i] = 0; 00371 } 00372 // if not inside, return distance to first hit 00373 gfxReal firstHit=-1.0; 00374 int firstOId = -1; 00375 if(globGeoInitDebug) errMsg("IIIstart"," isect "<<org<<" f"<<flags<<" acc"<<mAccurateGeoinit); 00376 00377 if(mAccurateGeoinit) { 00378 while(!done) { 00379 // find first inside intersection 00380 ntlTriangle *triIns = NULL; 00381 distance = -1.0; 00382 ntlVec3Gfx normal(0.0); 00383 if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true); 00384 else mpGiTree->intersect(ray,distance,normal, triIns, flags, true); 00385 if(triIns) { 00386 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance; 00387 LbmFloat orientation = dot(normal, dir); 00388 OId = triIns->getObjectId(); 00389 if(orientation<=0.0) { 00390 // outside hit 00391 normal *= -1.0; 00392 mGiObjInside[OId]++; 00393 if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = 1; 00394 if(globGeoInitDebug) errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation); 00395 } else { 00396 // inside hit 00397 mGiObjInside[OId]++; 00398 if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance; 00399 if(globGeoInitDebug) errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg<<" orient:"<<orientation); 00400 if(giObjFirstHistSide[OId]==0) giObjFirstHistSide[OId] = -1; 00401 } 00402 norg += normal * getVecEpsilon(); 00403 ray = ntlRay(norg, dir, 0, 1.0, mpGlob); 00404 // remember first hit distance, in case we're not 00405 // inside anything 00406 if(firstHit<0.0) { 00407 firstHit = distance; 00408 firstOId = OId; 00409 } 00410 } else { 00411 // no more intersections... return false 00412 done = true; 00413 } 00414 } 00415 00416 distance = -1.0; 00417 for(size_t i=0; i<mGiObjInside.size(); i++) { 00418 if(mGiObjInside[i]>0) { 00419 bool mess = false; 00420 if((mGiObjInside[i]%2)==1) { 00421 if(giObjFirstHistSide[i] != -1) mess=true; 00422 } else { 00423 if(giObjFirstHistSide[i] != 1) mess=true; 00424 } 00425 if(mess) { 00426 //errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] ); 00427 globGICPIProblems++; 00428 mGiObjInside[i]++; // believe first hit side... 00429 } 00430 } 00431 } 00432 for(size_t i=0; i<mGiObjInside.size(); i++) { 00433 if(globGeoInitDebug) errMsg("CHIII","i"<<i<<" ins="<<mGiObjInside[i]<<" t"<<mGiObjDistance[i]<<" d"<<distance); 00434 if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) { 00435 if( (distance<0.0) || // first intersection -> good 00436 ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one 00437 ) { 00438 distance = mGiObjDistance[i]; 00439 OId = i; 00440 inside = true; 00441 } 00442 } 00443 } 00444 if(!inside) { 00445 distance = firstHit; 00446 OId = firstOId; 00447 } 00448 if(globGeoInitDebug) errMsg("CHIII","i"<<inside<<" fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId); 00449 00450 return inside; 00451 } else { 00452 00453 // find first inside intersection 00454 ntlTriangle *triIns = NULL; 00455 distance = -1.0; 00456 ntlVec3Gfx normal(0.0); 00457 if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true); 00458 else mpGiTree->intersect(ray,distance,normal, triIns, flags, true); 00459 if(triIns) { 00460 // check outside intersect 00461 LbmFloat orientation = dot(normal, dir); 00462 if(orientation<=0.0) return false; 00463 00464 OId = triIns->getObjectId(); 00465 return true; 00466 } 00467 return false; 00468 } 00469 } 00470 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse) { 00471 // shift ve ctors to avoid rounding errors 00472 org += ntlVec3Gfx(0.0001); //? 00473 OId = -1; 00474 ntlRay ray(org, dir, 0, 1.0, mpGlob); 00475 //int insCnt = 0; 00476 bool done = false; 00477 bool inside = false; 00478 for(size_t i=0; i<mGiObjInside.size(); i++) { 00479 mGiObjInside[i] = 0; 00480 mGiObjDistance[i] = -1.0; 00481 mGiObjSecondDist[i] = -1.0; 00482 } 00483 // if not inside, return distance to first hit 00484 gfxReal firstHit=-1.0; 00485 int firstOId = -1; 00486 thinHit = false; 00487 00488 if(mAccurateGeoinit) { 00489 while(!done) { 00490 // find first inside intersection 00491 ntlTriangle *triIns = NULL; 00492 distance = -1.0; 00493 ntlVec3Gfx normal(0.0); 00494 mpGiTree->intersect(ray,distance,normal, triIns, flags, true); 00495 if(triIns) { 00496 ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance; 00497 LbmFloat orientation = dot(normal, dir); 00498 OId = triIns->getObjectId(); 00499 if(orientation<=0.0) { 00500 // outside hit 00501 normal *= -1.0; 00502 //mGiObjDistance[OId] = -1.0; 00503 //errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg); 00504 } else { 00505 // inside hit 00506 //if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance; 00507 //errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg); 00508 if(mGiObjInside[OId]==1) { 00509 // second inside hit 00510 if(mGiObjSecondDist[OId]<0.0) mGiObjSecondDist[OId] = distance; 00511 } 00512 } 00513 mGiObjInside[OId]++; 00514 // always store first hit for thin obj init 00515 if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance; 00516 00517 norg += normal * getVecEpsilon(); 00518 ray = ntlRay(norg, dir, 0, 1.0, mpGlob); 00519 // remember first hit distance, in case we're not 00520 // inside anything 00521 if(firstHit<0.0) { 00522 firstHit = distance; 00523 firstOId = OId; 00524 } 00525 } else { 00526 // no more intersections... return false 00527 done = true; 00528 //if(insCnt%2) inside=true; 00529 } 00530 } 00531 00532 distance = -1.0; 00533 // standard inside check 00534 for(size_t i=0; i<mGiObjInside.size(); i++) { 00535 if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) { 00536 if( (distance<0.0) || // first intersection -> good 00537 ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one 00538 ) { 00539 distance = mGiObjDistance[i]; 00540 OId = i; 00541 inside = true; 00542 } 00543 } 00544 } 00545 // now check for thin hits 00546 if(!inside) { 00547 distance = -1.0; 00548 for(size_t i=0; i<mGiObjInside.size(); i++) { 00549 if((mGiObjInside[i]>=2)&&(mGiObjDistance[i]>0.0)&&(mGiObjSecondDist[i]>0.0)&& 00550 (mGiObjDistance[i]<1.0*halfCellsize)&&(mGiObjSecondDist[i]<2.0*halfCellsize) ) { 00551 if( (distance<0.0) || // first intersection -> good 00552 ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one 00553 ) { 00554 distance = mGiObjDistance[i]; 00555 OId = i; 00556 inside = true; 00557 thinHit = true; 00558 } 00559 } 00560 } 00561 } 00562 if(!inside) { 00563 // check for hit in this cell, opposite to current dir (only recurse once) 00564 if(recurse) { 00565 gfxReal r_distance; 00566 int r_OId; 00567 bool ret = geoInitCheckPointInside(org, dir*-1.0, flags, r_OId, r_distance, halfCellsize, thinHit, false); 00568 if((ret)&&(thinHit)) { 00569 OId = r_OId; 00570 distance = 0.0; 00571 return true; 00572 } 00573 } 00574 } 00575 // really no hit... 00576 if(!inside) { 00577 distance = firstHit; 00578 OId = firstOId; 00579 /*if((mGiObjDistance[OId]>0.0)&&(mGiObjSecondDist[OId]>0.0)) { 00580 const gfxReal thisdist = mGiObjSecondDist[OId]-mGiObjDistance[OId]; 00581 // dont walk over this cell... 00582 if(thisdist<halfCellsize) distance-=2.0*halfCellsize; 00583 } // ? */ 00584 } 00585 //errMsg("CHIII","i"<<inside<<" fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId); 00586 00587 return inside; 00588 } else { 00589 00590 // find first inside intersection 00591 ntlTriangle *triIns = NULL; 00592 distance = -1.0; 00593 ntlVec3Gfx normal(0.0); 00594 mpGiTree->intersect(ray,distance,normal, triIns, flags, true); 00595 if(triIns) { 00596 // check outside intersect 00597 LbmFloat orientation = dot(normal, dir); 00598 if(orientation<=0.0) return false; 00599 00600 OId = triIns->getObjectId(); 00601 return true; 00602 } 00603 return false; 00604 } 00605 } 00606 00607 00608 /*****************************************************************************/ 00609 ntlVec3Gfx LbmSolverInterface::getGeoMaxMovementVelocity(LbmFloat simtime, LbmFloat stepsize) { 00610 ntlVec3Gfx max(0.0); 00611 if(mpGlob == NULL) return max; 00612 // mpGiObjects has to be inited here... 00613 00614 for(int i=0; i< (int)mpGiObjects->size(); i++) { 00615 //errMsg("LbmSolverInterface::getGeoMaxMovementVelocity","i="<<i<<" "<< (*mpGiObjects)[i]->getName() ); // DEBUG 00616 if( (*mpGiObjects)[i]->getGeoInitType() & (FGI_FLUID|FGI_MBNDINFLOW) ){ 00617 //ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime); 00618 ntlVec3Gfx orgvel = (*mpGiObjects)[i]->calculateMaxVel( simtime, simtime+stepsize ); 00619 if( normNoSqrt(orgvel) > normNoSqrt(max) ) { max = orgvel; } 00620 //errMsg("MVT","i"<<i<<" a"<< (*mpGiObjects)[i]->getInitialVelocity(simtime)<<" o"<<orgvel ); // DEBUG 00621 // TODO check if inflow and simtime 00622 ntlVec3Gfx inivel = (*mpGiObjects)[i]->getInitialVelocity(simtime); 00623 if( normNoSqrt(inivel) > normNoSqrt(max) ) { max = inivel; } 00624 } 00625 } 00626 errMsg("LbmSolverInterface::getGeoMaxMovementVelocity", "max="<< max ); // DEBUG 00627 return max; 00628 } 00629 00630 00631 00632 00633 /*******************************************************************************/ 00635 /*******************************************************************************/ 00636 00637 00638 00639 00640 /*****************************************************************************/ 00642 void 00643 LbmSolverInterface::addCellToMarkedList( CellIdentifierInterface *cid ) { 00644 for(size_t i=0; i<mMarkedCells.size(); i++) { 00645 // check if cids alreay in 00646 if( mMarkedCells[i]->equal(cid) ) return; 00647 //mMarkedCells[i]->setEnd(false); 00648 } 00649 mMarkedCells.push_back( cid ); 00650 //cid->setEnd(true); 00651 } 00652 00653 /*****************************************************************************/ 00655 CellIdentifierInterface* 00656 LbmSolverInterface::markedGetFirstCell( ) { 00657 if(mMarkedCells.size() > 0){ return mMarkedCells[0]; } 00658 return NULL; 00659 } 00660 00661 CellIdentifierInterface* 00662 LbmSolverInterface::markedAdvanceCell() { 00663 mMarkedCellIndex++; 00664 if(mMarkedCellIndex>=(int)mMarkedCells.size()) return NULL; 00665 return mMarkedCells[mMarkedCellIndex]; 00666 } 00667 00668 void LbmSolverInterface::markedClearList() { 00669 // FIXME free cids? 00670 mMarkedCells.clear(); 00671 } 00672 00673 /*******************************************************************************/ 00675 /*******************************************************************************/ 00676 00677 00678 00679 // 32k 00680 string convertSingleFlag2String(CellFlagType cflag) { 00681 CellFlagType flag = cflag; 00682 if(flag == CFUnused ) return string("cCFUnused"); 00683 if(flag == CFEmpty ) return string("cCFEmpty"); 00684 if(flag == CFBnd ) return string("cCFBnd"); 00685 if(flag == CFBndNoslip ) return string("cCFBndNoSlip"); 00686 if(flag == CFBndFreeslip ) return string("cCFBndFreeSlip"); 00687 if(flag == CFBndPartslip ) return string("cCFBndPartSlip"); 00688 if(flag == CFNoInterpolSrc ) return string("cCFNoInterpolSrc"); 00689 if(flag == CFFluid ) return string("cCFFluid"); 00690 if(flag == CFInter ) return string("cCFInter"); 00691 if(flag == CFNoNbFluid ) return string("cCFNoNbFluid"); 00692 if(flag == CFNoNbEmpty ) return string("cCFNoNbEmpty"); 00693 if(flag == CFNoDelete ) return string("cCFNoDelete"); 00694 if(flag == CFNoBndFluid ) return string("cCFNoBndFluid"); 00695 if(flag == CFGrNorm ) return string("cCFGrNorm"); 00696 if(flag == CFGrFromFine ) return string("cCFGrFromFine"); 00697 if(flag == CFGrFromCoarse ) return string("cCFGrFromCoarse"); 00698 if(flag == CFGrCoarseInited ) return string("cCFGrCoarseInited"); 00699 if(flag == CFMbndInflow ) return string("cCFMbndInflow"); 00700 if(flag == CFMbndOutflow ) return string("cCFMbndOutflow"); 00701 if(flag == CFInvalid ) return string("cfINVALID"); 00702 00703 std::ostringstream mult; 00704 int val = 0; 00705 if(flag != 0) { 00706 while(! (flag&1) ) { 00707 flag = flag>>1; 00708 val++; 00709 } 00710 } else { 00711 val = -1; 00712 } 00713 if(val>=24) { 00714 mult << "cfOID_" << (flag>>24) <<"_TYPE"; 00715 } else { 00716 mult << "cfUNKNOWN_" << val <<"_TYPE"; 00717 } 00718 return mult.str(); 00719 } 00720 00722 string convertCellFlagType2String( CellFlagType cflag ) { 00723 int flag = cflag; 00724 00725 const int jmax = sizeof(CellFlagType)*8; 00726 bool somefound = false; 00727 std::ostringstream mult; 00728 mult << "["; 00729 for(int j=0; j<jmax ; j++) { 00730 if(flag& (1<<j)) { 00731 if(somefound) mult << "|"; 00732 mult << j<<"<"<< convertSingleFlag2String( (CellFlagType)(1<<j) ); // this call should always be _non_-recursive 00733 somefound = true; 00734 } 00735 }; 00736 mult << "]"; 00737 00738 // return concatenated string 00739 if(somefound) return mult.str(); 00740 00741 // empty? 00742 return string("[emptyCFT]"); 00743 } 00744