diff options
author | Nils Thuerey <nils@thuerey.de> | 2006-06-12 16:55:51 +0400 |
---|---|---|
committer | Nils Thuerey <nils@thuerey.de> | 2006-06-12 16:55:51 +0400 |
commit | b0a22902ec1b6586325a027358415f5bdcb828f9 (patch) | |
tree | cc557b32778bf7406f112447a0a54c8de40a6aed /intern/elbeem | |
parent | a0d94e672730e854d05157133a755ee3257f373a (diff) |
- another minor solver update to fix
obstacle fluid surface generation bug
- also contains some code clean ups
and safer initializations
Diffstat (limited to 'intern/elbeem')
-rw-r--r-- | intern/elbeem/intern/attributes.cpp | 37 | ||||
-rw-r--r-- | intern/elbeem/intern/elbeem.cpp | 21 | ||||
-rw-r--r-- | intern/elbeem/intern/ntl_geometrymodel.cpp | 46 | ||||
-rw-r--r-- | intern/elbeem/intern/ntl_geometrymodel.h | 13 | ||||
-rw-r--r-- | intern/elbeem/intern/ntl_ray.h | 2 | ||||
-rw-r--r-- | intern/elbeem/intern/particletracer.cpp | 6 | ||||
-rw-r--r-- | intern/elbeem/intern/simulation_object.cpp | 4 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_class.h | 18 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_init.cpp | 14 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_interface.cpp | 7 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_main.cpp | 44 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_relax.h | 8 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_util.cpp | 86 |
13 files changed, 205 insertions, 101 deletions
diff --git a/intern/elbeem/intern/attributes.cpp b/intern/elbeem/intern/attributes.cpp index f970f6ee5af..b7f194582d1 100644 --- a/intern/elbeem/intern/attributes.cpp +++ b/intern/elbeem/intern/attributes.cpp @@ -472,35 +472,45 @@ bool AttributeList::ignoreParameter(string name, string source) { // read channels AnimChannel<int> AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel<int>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<int>(defaultValue); } AnimChannel<int> ret = find(name)->getChannelInt(); find(name)->setUsed(true); channelSimplifyi(ret); return ret; } AnimChannel<double> AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel<double>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<double>(defaultValue); } AnimChannel<double> ret = find(name)->getChannelFloat(); find(name)->setUsed(true); channelSimplifyd(ret); return ret; } AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel<ntlVec3d>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<ntlVec3d>(defaultValue); } AnimChannel<ntlVec3d> ret = find(name)->getChannelVec3d(); find(name)->setUsed(true); channelSimplifyVd(ret); return ret; } AnimChannel<ntlSetVec3f> AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel<ntlSetVec3f>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<ntlSetVec3f>(defaultValue); } AnimChannel<ntlSetVec3f> ret = find(name)->getChannelSetVec3f(); find(name)->setUsed(true); //channelSimplifyVf(ret); return ret; } AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel<float>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<float>(defaultValue); } AnimChannel<double> convert = find(name)->getChannelFloat(); find(name)->setUsed(true); channelSimplifyd(convert); @@ -514,7 +524,9 @@ AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float return ret; } AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel<ntlVec3f>(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); } + return AnimChannel<ntlVec3f>(defaultValue); } AnimChannel<ntlVec3d> convert = find(name)->getChannelVec3d(); // convert to float @@ -737,6 +749,19 @@ void AnimChannel<Scalar>::debugPrintChannel() { } +// hack to force instantiation +void __forceAnimChannelInstantiation() { + AnimChannel< float > tmp1; + AnimChannel< double > tmp2; + AnimChannel< string > tmp3; + AnimChannel< ntlVector3Dim<float> > tmp4; + tmp1.debugPrintChannel(); + tmp2.debugPrintChannel(); + tmp3.debugPrintChannel(); + tmp4.debugPrintChannel(); +} + + ntlSetVec3f::ntlSetVec3f(double v ) { mVerts.clear(); mVerts.push_back( ntlVec3f(v) ); diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp index 907401395c5..d4242da4d5e 100644 --- a/intern/elbeem/intern/elbeem.cpp +++ b/intern/elbeem/intern/elbeem.cpp @@ -115,32 +115,45 @@ void elbeemGetErrorString(char *buffer) { extern "C" void elbeemResetMesh(elbeemMesh *mesh) { if(!mesh) return; + // init typedef struct elbeemMesh mesh->type = 0; - mesh->parentDomainId = 0; + + mesh->parentDomainId = 0; + + /* vertices */ mesh->numVertices = 0; mesh->vertices = NULL; + + mesh->channelSizeVertices = 0; + mesh->channelVertices = NULL; + + /* triangles */ mesh->numTriangles = 0; mesh->triangles = NULL; + /* animation channels */ mesh->channelSizeTranslation = 0; mesh->channelTranslation = NULL; mesh->channelSizeRotation = 0; mesh->channelRotation = NULL; mesh->channelSizeScale = 0; mesh->channelScale = NULL; + + /* active channel */ mesh->channelSizeActive = 0; mesh->channelActive = NULL; + mesh->channelSizeInitialVel = 0; mesh->channelInitialVel = NULL; + mesh->localInivelCoords = 0; mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP; - mesh->volumeInitType= OB_VOLUMEINIT_VOLUME; mesh->obstaclePartslip= 0.; - mesh->channelSizeVertices = 0; - mesh->channelVertices = NULL; + mesh->volumeInitType= OB_VOLUMEINIT_VOLUME; + /* name of the mesh, mostly for debugging */ mesh->name = "[unnamed]"; } diff --git a/intern/elbeem/intern/ntl_geometrymodel.cpp b/intern/elbeem/intern/ntl_geometrymodel.cpp index 16d18c4b238..ba6e638f507 100644 --- a/intern/elbeem/intern/ntl_geometrymodel.cpp +++ b/intern/elbeem/intern/ntl_geometrymodel.cpp @@ -28,9 +28,7 @@ ntlGeometryObjModel::ntlGeometryObjModel( void ) : mLoaded( false ), mTriangles(), mVertices(), mNormals(), mcAniVerts(), mcAniNorms(), - mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.), - mvCPSStart(-10000.), mvCPSEnd(10000.), - mCPSFilename(""), mCPSWidth(0.1), mCPSTimestep(1.) + mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.) { } @@ -52,6 +50,39 @@ bool ntlGeometryObjModel::getMeshAnimated() { return ret; } +/*! calculate max extends of (ani) mesh */ +void ntlGeometryObjModel::getExtends(ntlVec3Gfx &sstart, ntlVec3Gfx &send) { + bool ini=false; + ntlVec3Gfx start(0.),end(0.); + for(int s=0; s<=(int)mcAniVerts.accessValues().size(); s++) { + vector<ntlVec3f> *sverts; + if(mcAniVerts.accessValues().size()>0) { + if(s==(int)mcAniVerts.accessValues().size()) continue; + sverts = &(mcAniVerts.accessValues()[s].mVerts); + } else sverts = &mVertices; + + for(int i=0; i<(int)sverts->size(); i++) { + + if(!ini) { + start=(*sverts)[i]; + end=(*sverts)[i]; + //errMsg("getExtends","ini "<<s<<","<<i<<" "<<start<<","<<end); + ini=true; + } else { + for(int j=0; j<3; j++) { + if(start[j] > (*sverts)[i][j]) { start[j]= (*sverts)[i][j]; } + if(end[j] < (*sverts)[i][j]) { end[j] = (*sverts)[i][j]; } + } + //errMsg("getExtends","check "<<s<<","<<i<<" "<<start<<","<<end<<" "<< (*sverts)[i]); + } + + } + } + sstart=start; + send=end; +} + + /*****************************************************************************/ /* Init attributes etc. of this object */ /*****************************************************************************/ @@ -90,11 +121,6 @@ void ntlGeometryObjModel::initialize(ntlRenderGlobals *glob) mAniTimeScale = mpAttrs->readFloat("ani_timescale", mAniTimeScale,"ntlGeometryObjModel", "mAniTimeScale", false); mAniTimeOffset = mpAttrs->readFloat("ani_timeoffset", mAniTimeOffset,"ntlGeometryObjModel", "mAniTimeOffset", false); - mCPSWidth = mpAttrs->readFloat("cps_width", mCPSWidth,"ntlGeometryObjModel", "mCPSWidth", false); - mCPSTimestep = mpAttrs->readFloat("cps_timestep", mCPSTimestep,"ntlGeometryObjModel", "mCPSTimestep", false); - mvCPSStart = vec2G( mpAttrs->readVec3d("cps_start", vec2D(mvCPSStart),"ntlGeometryObjModel", "mvCPSStart", false)); - mvCPSEnd = vec2G( mpAttrs->readVec3d("cps_end", vec2D(mvCPSEnd),"ntlGeometryObjModel", "mvCPSEnd", false)); - // continue with standard obj if(loadBobjModel(mFilename)==0) mLoaded=1; if(!mLoaded) { @@ -322,7 +348,7 @@ int ntlGeometryObjModel::loadBobjModel(string filename) bytesRead += gzread(gzf, &frameTime, sizeof(frameTime) ); //if(bytesRead!=3*sizeof(float)) { if(bytesRead!=sizeof(float)) { - debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' no ani sets. ", 10 ); + debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' end of gzfile. ", 10 ); if(anitimes.size()>0) { // finally init channels and stop reading file mcAniVerts = AnimChannel<ntlSetVec3f>(aniverts,anitimes); @@ -428,7 +454,7 @@ ntlGeometryObjModel::getTriangles(double t, vector<ntlTriangle> *triangles, (*normals)[startvert+i] = mNormals[i]; } - triangles->reserve(triangles->size() + mTriangles.size() ); + triangles->reserve(triangles->size() + mTriangles.size()/3 ); for(int i=0; i<(int)mTriangles.size(); i+=3) { int trip[3]; trip[0] = startvert+mTriangles[i+0]; diff --git a/intern/elbeem/intern/ntl_geometrymodel.h b/intern/elbeem/intern/ntl_geometrymodel.h index f911279489d..95cadc47925 100644 --- a/intern/elbeem/intern/ntl_geometrymodel.h +++ b/intern/elbeem/intern/ntl_geometrymodel.h @@ -44,8 +44,8 @@ class ntlGeometryObjModel : public ntlGeometryObject /*! init triangle divisions */ virtual void calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri); - /*! do ani mesh CPS */ - void calculateCPS(string filename); + /*! calculate max extends of (ani) mesh */ + void getExtends(ntlVec3Gfx &start, ntlVec3Gfx &end); private: @@ -71,12 +71,6 @@ class ntlGeometryObjModel : public ntlGeometryObject /*! timing mapping & offset for config files */ double mAniTimeScale, mAniTimeOffset; - /*! ani mesh cps params */ - ntlVec3Gfx mvCPSStart, mvCPSEnd; - string mCPSFilename; - gfxReal mCPSWidth, mCPSTimestep; - - public: /* Access methods */ @@ -87,6 +81,9 @@ class ntlGeometryObjModel : public ntlGeometryObject inline ntlVec3Gfx getEnd( void ){ return mvEnd; } inline void setEnd( const ntlVec3Gfx &set ){ mvEnd = set; } + inline bool getLoaded( void ){ return mLoaded; } + inline void setLoaded( bool set ){ mLoaded = set; } + /*! set data file name */ inline void setFilename(string set) { mFilename = set; } }; diff --git a/intern/elbeem/intern/ntl_ray.h b/intern/elbeem/intern/ntl_ray.h index 2321cee15e9..8c0188aee30 100644 --- a/intern/elbeem/intern/ntl_ray.h +++ b/intern/elbeem/intern/ntl_ray.h @@ -331,6 +331,8 @@ public: mGeos.push_back( geo ); geo->setObjectId(mGeos.size()); } + /*! Add a geo object to the scene, warning - only needed for hand init */ + inline void addGeoObject(ntlGeometryObject *geo) { mObjects.push_back( geo ); } /*! Acces a certain object */ inline ntlGeometryObject *getObject(int id) { diff --git a/intern/elbeem/intern/particletracer.cpp b/intern/elbeem/intern/particletracer.cpp index c49c1b1a07e..1d15cecfd66 100644 --- a/intern/elbeem/intern/particletracer.cpp +++ b/intern/elbeem/intern/particletracer.cpp @@ -154,7 +154,7 @@ void ParticleTracer::cleanup() { *! dump particles if desired *****************************************************************************/ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) { - debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts, 10); // DEBUG + debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG if( (dumptype==DUMP_FULLGEOMETRY)&& @@ -305,7 +305,7 @@ void ParticleTracer::checkTrails(double time) { /****************************************************************************** * Get triangles for rendering *****************************************************************************/ -void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles, +void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, vector<ntlVec3Gfx> *vertices, vector<ntlVec3Gfx> *normals, int objectId ) { @@ -326,7 +326,7 @@ void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles, int segments = mPartSegments; ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2])); ntlVec3Gfx trans = mStart; - t = 0.; // doesnt matter + time = 0.; // doesnt matter for(size_t t=0; t<mPrevs.size()+1; t++) { vector<ParticleObject> *dparts; diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp index 5d8a73e6ab8..e34030180a5 100644 --- a/intern/elbeem/intern/simulation_object.cpp +++ b/intern/elbeem/intern/simulation_object.cpp @@ -96,6 +96,7 @@ void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) { *mpElbeemSettings = *settings; mGeoInitId = settings->domainId+1; + debMsgStd("SimulationObject",DM_MSG,"mGeoInitId="<<mGeoInitId<<", domainId="<<settings->domainId, 8); } /****************************************************************************** @@ -115,7 +116,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) } - this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"LbmSolverInterface", "mGeoInitId", false); + mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"LbmSolverInterface", "mGeoInitId", false); //mDimension, mSolverType are deprecated string mSolverType(""); mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false ); @@ -297,6 +298,7 @@ void SimulationObject::step( void ) // dont advance for stopped time mpLbm->step(); mTime += mpParam->getTimestep(); +//if(mTime>0.001) { errMsg("DEBUG!!!!!!!!","quit mlsu..."); exit(1); } // PROFILE DEBUG TEST! } if(mpLbm->getPanic()) mPanic = true; diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h index 10b172dde4b..46462ec52de 100644 --- a/intern/elbeem/intern/solver_class.h +++ b/intern/elbeem/intern/solver_class.h @@ -106,10 +106,6 @@ #endif #endif -#if LBM_INCLUDE_TESTSOLVERS==1 -#include "solver_test.h" -#endif // LBM_INCLUDE_TESTSOLVERS==1 - /*****************************************************************************/ /*! cell access classes */ class UniformFsgrCellIdentifier : @@ -467,20 +463,6 @@ class LbmFsgrSolver : # endif // FSGR_STRICT_DEBUG==1 bool mUseTestdata; -#if LBM_INCLUDE_TESTSOLVERS==1 - // test functions - LbmTestdata *mpTest; - void initTestdata(); - void destroyTestdata(); - void handleTestdata(); - void set3dHeight(int ,int ); - - void initCpdata(); - void handleCpdata(); - public: - // needed for testdata - void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz); -#endif // LBM_INCLUDE_TESTSOLVERS==1 public: // former LbmModelLBGK functions // relaxation funtions - implemented together with relax macros diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp index 908355fdb13..226af4ec76d 100644 --- a/intern/elbeem/intern/solver_init.cpp +++ b/intern/elbeem/intern/solver_init.cpp @@ -411,7 +411,7 @@ LbmFsgrSolver::~LbmFsgrSolver() { if(!this->mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; } #if COMPRESSGRIDS==1 - delete mLevel[mMaxRefine].mprsCells[1]; + delete [] mLevel[mMaxRefine].mprsCells[1]; mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL; #endif // COMPRESSGRIDS==1 @@ -469,6 +469,12 @@ void LbmFsgrSolver::parseAttrList() // FIXME check needed? mFVArea = this->mpAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false ); + // debugging - skip some time... + double starttimeskip = 0.; + starttimeskip = this->mpAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false ); + mSimulationTime += starttimeskip; + if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1); + #if LBM_INCLUDE_TESTSOLVERS==1 mUseTestdata = 0; mUseTestdata = this->mpAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false); @@ -483,7 +489,7 @@ void LbmFsgrSolver::parseAttrList() mUseTestdata = 0; if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check #endif // LBM_INCLUDE_TESTSOLVERS!=1 - if(mUseTestdata) { mMaxRefine=0; } // force fsgr off + if(mUseTestdata) { mMaxRefine=0; } // force fsgr off if(mMaxRefine==0) mInitialCsmago=0.02; mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false ); @@ -1184,6 +1190,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip; } break; + // off */ /* case FGI_BNDPART: rhomass = BND_FILL; otype = ntype = CFBnd|CFBndPartslip; @@ -1191,7 +1198,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { case FGI_BNDFREE: rhomass = BND_FILL; otype = ntype = CFBnd|CFBndFreeslip; break; - // */ + // off */ case FGI_BNDNO: rhomass = BND_FILL; otype = ntype = CFBnd|CFBndNoslip; break; @@ -1876,7 +1883,6 @@ void LbmFsgrSolver::initFreeSurfaces() { mLevel[lev].setCurr ^= 1; } // copy back...? - } /*****************************************************************************/ diff --git a/intern/elbeem/intern/solver_interface.cpp b/intern/elbeem/intern/solver_interface.cpp index d682445645d..c4e8d22a173 100644 --- a/intern/elbeem/intern/solver_interface.cpp +++ b/intern/elbeem/intern/solver_interface.cpp @@ -44,6 +44,7 @@ LbmSolverInterface::LbmSolverInterface() : mName("lbm_default") , mpIso( NULL ), mIsoValue(0.499), mSilent(false) , + mLbmInitId(1) , mpGiTree( NULL ), mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ), mRefinementDesired(0), @@ -284,7 +285,7 @@ void LbmSolverInterface::initGeoTree() { if(mpGiTree != NULL) delete mpGiTree; char treeFlag = (1<<(this->mLbmInitId+4)); mpGiTree = new ntlTree( - 15, 8, // warning - fixed values for depth & maxtriangles here... + 15, 8, // TREEwarning - fixed values for depth & maxtriangles here... scene, treeFlag ); } @@ -299,6 +300,7 @@ void LbmSolverInterface::freeGeoTree() { int globGeoInitDebug = 0; +int globGICPIProblems = 0; /*****************************************************************************/ /*! check for a certain flag type at position org */ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance) { @@ -371,7 +373,8 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int if(giObjFirstHistSide[i] != 1) mess=true; } if(mess) { - errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] ); + //errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] ); + globGICPIProblems++; mGiObjInside[i]++; // believe first hit side... } } diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp index 2888c7e172e..8fa3b8d8a65 100644 --- a/intern/elbeem/intern/solver_main.cpp +++ b/intern/elbeem/intern/solver_main.cpp @@ -294,6 +294,8 @@ LbmFsgrSolver::mainLoop(int lev) int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine); //{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIMcDimension<<" mlsz:"<< mLevel[mMaxRefine].lSizez<<" zmax1:"<<getForZMax1(mMaxRefine) ); } // DEBUG #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz); +#define LIST_EMPTY(x) mListEmpty.push_back( x ); +#define LIST_FULL(x) mListFull.push_back( x ); #endif // PARALLEL==1 // local to loop @@ -360,22 +362,15 @@ LbmFsgrSolver::mainLoop(int lev) } // COMPRT #if PARALLEL==0 - const int id = 0, Nthrds = 1; -#endif // PARALLEL==1 - const int Nj = mLevel[mMaxRefine].lSizey; - int jstart = 0+( id * (Nj / Nthrds) ); - int jend = 0+( (id+1) * (Nj / Nthrds) ); - if( ((Nj/Nthrds) *Nthrds) != Nj) { - errMsg("LbmFsgrSolver","Invalid domain size Nj="<<Nj<<" Nthrds="<<Nthrds); - } - // cutoff obstacle boundary - if(jstart<1) jstart = 1; - if(jend>mLevel[mMaxRefine].lSizey-1) jend = mLevel[mMaxRefine].lSizey-1; - -#if PARALLEL==1 + //const int id = 0, Nthrds = 1; + const int jstart = 0; + const int jend = mLevel[mMaxRefine].lSizey; +//#endif // PARALLEL==1 +#else // PARALLEL==1 PARA_INITIALIZE(); - //errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug + errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug #endif // PARALLEL==1 + for(int k=kstart;k!=kend;k+=kdir) { //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug @@ -459,13 +454,14 @@ LbmFsgrSolver::mainLoop(int lev) changeFlag(lev, i,j,k, TSET(lev), CFInter); // same as ifemptied for if below - LbmPoint emptyp; emptyp.flag = 0; - emptyp.x = i; emptyp.y = j; emptyp.z = k; + LbmPoint oemptyp; oemptyp.flag = 0; + oemptyp.x = i; oemptyp.y = j; oemptyp.z = k; #if PARALLEL==1 - calcListEmpty[id].push_back( emptyp ); + //calcListEmpty[id].push_back( oemptyp ); #else // PARALLEL==1 - mListEmpty.push_back( emptyp ); + //mListEmpty.push_back( oemptyp ); #endif // PARALLEL==1 + LIST_EMPTY(oemptyp); calcCellsEmptied++; continue; } @@ -601,6 +597,7 @@ LbmFsgrSolver::mainLoop(int lev) // for fluid cells - just the f_i difference from streaming to empty cells ---- numRecons = 0; bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip)); + //onlyBndnb = false; // DEBUG test off FORDF1 { // dfl loop recons[l] = 0; @@ -1008,6 +1005,7 @@ LbmFsgrSolver::mainLoop(int lev) // looks much nicer... LISTTRICK #if FSGR_LISTTRICK==1 + if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; } if(newFlag&CFNoBndFluid) { // test NEW TEST if(!iffilled) { // remove cells independent from amount of change... @@ -1035,10 +1033,11 @@ LbmFsgrSolver::mainLoop(int lev) if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1; // NEWSURFT filledp.x = i; filledp.y = j; filledp.z = k; #if PARALLEL==1 - calcListFull[id].push_back( filledp ); + //calcListFull[id].push_back( filledp ); #else // PARALLEL==1 - mListFull.push_back( filledp ); + //mListFull.push_back( filledp ); #endif // PARALLEL==1 + LIST_FULL(filledp); //this->mNumFilledCells++; // DEBUG calcCellsFilled++; } @@ -1047,10 +1046,11 @@ LbmFsgrSolver::mainLoop(int lev) if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; // NEWSURFT emptyp.x = i; emptyp.y = j; emptyp.z = k; #if PARALLEL==1 - calcListEmpty[id].push_back( emptyp ); + //calcListEmpty[id].push_back( emptyp ); #else // PARALLEL==1 - mListEmpty.push_back( emptyp ); + //mListEmpty.push_back( emptyp ); #endif // PARALLEL==1 + LIST_EMPTY(emptyp); //this->mNumEmptiedCells++; // DEBUG calcCellsEmptied++; } diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h index eba6742169e..1a1e2c01172 100644 --- a/intern/elbeem/intern/solver_relax.h +++ b/intern/elbeem/intern/solver_relax.h @@ -17,7 +17,6 @@ -#if LBM_INCLUDE_TESTSOLVERS!=1 // off for non testing #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \ @@ -25,11 +24,8 @@ uy += (grav)[1]; \ uz += (grav)[2]; \ -#else // LBM_INCLUDE_TESTSOLVERS!=1 +#define TEST_IF_CHECK -// defined in test.h - -#endif // LBM_INCLUDE_TESTSOLVERS!=1 /****************************************************************************** @@ -1125,7 +1121,7 @@ inline void LbmFsgrSolver::collideArrays( for(l=0; l<this->cDfNum; l++) { df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; } - if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df); + //if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df); mux = ux; muy = uy; diff --git a/intern/elbeem/intern/solver_util.cpp b/intern/elbeem/intern/solver_util.cpp index 8346299b39a..05cfcaa8e8e 100644 --- a/intern/elbeem/intern/solver_util.cpp +++ b/intern/elbeem/intern/solver_util.cpp @@ -53,33 +53,60 @@ void LbmFsgrSolver::prepareVisualization( void ) { for(int j=1;j<mLevel[lev].lSizey-1;j++) for(int i=1;i<mLevel[lev].lSizex-1;i++) { const CellFlagType cflag = RFLAG(lev, i,j,k,workSet); - //continue; // OFF DEBUG - if(cflag&(CFBnd|CFEmpty)) { + //if(cflag&(CFBnd|CFEmpty)) { + if(cflag&(CFBnd)) { continue; - } else if( (cflag&CFInter) ) { - //} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) && (cflag&CFNoNbFluid) ) { - //} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) ) { + } else if( (cflag&CFEmpty) ) { + //continue; // OFF DEBUG + int noslipbnd = 0; + int intercnt = 0; + CellFlagType nbored; + FORDF1 { + const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l); + if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; } + if(nbflag&CFInter){ intercnt++; } + nbored |= nbflag; + } + //? val = (QCELL(lev, i,j,k,workSet, dFfrac)); + if((noslipbnd)&&(intercnt>6)) { + //if(val<minval) val = minval; + //*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); + *this->mpIso->lbmGetData(i,j,ZKOFF) += minval; + } else if((noslipbnd)&&(intercnt>0)) { + *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.95; + } else { + } + continue; + + //} else if( (cflag&CFInter) ) { + + } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) { + // optimized fluid + val = 1.; + + } else if( (cflag&(CFInter|CFFluid)) ) { int noslipbnd = 0; FORDF1 { const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l); if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(CFBndNoslip)){ noslipbnd=1; l=100; } } + // no empty nb interface cells are treated as full + if(cflag&(CFNoNbEmpty|CFFluid)) { + val=1.0; + } val = (QCELL(lev, i,j,k,workSet, dFfrac)); + if(noslipbnd) { //errMsg("NEWVAL", PRINT_IJK<<" val"<<val <<" f"<< convertCellFlagType2String(cflag)<<" "<<noslipbnd); //" nbfl"<<convertCellFlagType2String(nbored) ); if(val<minval) val = minval; - *this->mpIso->lbmGetData( i , j ,ZKOFF ) += minval-( val * mIsoWeight[13] ); + *this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); } // */ - - // no empty nb interface cells are treated as full - if(cflag&CFNoNbEmpty) { - val=1.0; - } - } else { // fluid? - val = 1.0; + } else { // unused? + continue; + // old fluid val = 1.0; } // */ *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); @@ -885,10 +912,11 @@ void LbmFsgrSolver::advanceParticles() { } void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) { + int workSet = mLevel[mMaxRefine].setCurr; + std::ostringstream name; + // debug - raw dump of ffrac values if(getenv("ELBEEM_RAWDEBUGDUMP")) { - int workSet = mLevel[mMaxRefine].setCurr; - std::ostringstream name; //name <<"fill_" << this->mStepCnt <<".dump"; name << outfilename<< frameNrStr <<".dump"; FILE *file = fopen(name.str().c_str(),"w"); @@ -898,9 +926,12 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) { for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) { float val = 0.; - if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); + if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) { + val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); + if(val<0.) val=0.; + if(val>1.) val=1.; + } if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; - //fwrite( &val, sizeof(val), 1, file); // binary fprintf(file, "%f ",val); // text //errMsg("W", PRINT_IJK<<" val:"<<val); } @@ -912,6 +943,27 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt } // file } // */ + if(getenv("ELBEEM_BINDEBUGDUMP")) { + name << outfilename<< frameNrStr <<".bdump"; + FILE *file = fopen(name.str().c_str(),"w"); + if(file) { + for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { + for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) { + for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) { + float val = 0.; + if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) { + val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); + if(val<0.) val=0.; + if(val>1.) val=1.; + } + if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; + fwrite( &val, sizeof(val), 1, file); // binary + } + } + } + fclose(file); + } // file + } dumptype = 0; frameNr = 0; // get rid of warning } |