From b0a22902ec1b6586325a027358415f5bdcb828f9 Mon Sep 17 00:00:00 2001 From: Nils Thuerey Date: Mon, 12 Jun 2006 12:55:51 +0000 Subject: - another minor solver update to fix obstacle fluid surface generation bug - also contains some code clean ups and safer initializations --- intern/elbeem/intern/attributes.cpp | 37 ++++++++++--- intern/elbeem/intern/elbeem.cpp | 21 ++++++-- intern/elbeem/intern/ntl_geometrymodel.cpp | 46 ++++++++++++---- intern/elbeem/intern/ntl_geometrymodel.h | 13 ++--- intern/elbeem/intern/ntl_ray.h | 2 + intern/elbeem/intern/particletracer.cpp | 6 +-- intern/elbeem/intern/simulation_object.cpp | 4 +- intern/elbeem/intern/solver_class.h | 18 ------- intern/elbeem/intern/solver_init.cpp | 14 +++-- intern/elbeem/intern/solver_interface.cpp | 7 ++- intern/elbeem/intern/solver_main.cpp | 44 +++++++-------- intern/elbeem/intern/solver_relax.h | 8 +-- 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 AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<(defaultValue); } AnimChannel ret = find(name)->getChannelInt(); find(name)->setUsed(true); channelSimplifyi(ret); return ret; } AnimChannel AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<(defaultValue); } AnimChannel ret = find(name)->getChannelFloat(); find(name)->setUsed(true); channelSimplifyd(ret); return ret; } AnimChannel AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<(defaultValue); } AnimChannel ret = find(name)->getChannelVec3d(); find(name)->setUsed(true); channelSimplifyVd(ret); return ret; } AnimChannel AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<(defaultValue); } AnimChannel ret = find(name)->getChannelSetVec3f(); find(name)->setUsed(true); //channelSimplifyVf(ret); return ret; } AnimChannel AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<(defaultValue); } AnimChannel convert = find(name)->getChannelFloat(); find(name)->setUsed(true); channelSimplifyd(convert); @@ -514,7 +524,9 @@ AnimChannel AttributeList::readChannelSinglePrecFloat(string name, float return ret; } AnimChannel AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { return AnimChannel(defaultValue); } + if(!exists(name)) { + if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<(defaultValue); } AnimChannel convert = find(name)->getChannelVec3d(); // convert to float @@ -737,6 +749,19 @@ void AnimChannel::debugPrintChannel() { } +// hack to force instantiation +void __forceAnimChannelInstantiation() { + AnimChannel< float > tmp1; + AnimChannel< double > tmp2; + AnimChannel< string > tmp3; + AnimChannel< ntlVector3Dim > 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 *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 "< (*sverts)[i][j]) { start[j]= (*sverts)[i][j]; } + if(end[j] < (*sverts)[i][j]) { end[j] = (*sverts)[i][j]; } + } + //errMsg("getExtends","check "<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 '"<0) { // finally init channels and stop reading file mcAniVerts = AnimChannel(aniverts,anitimes); @@ -428,7 +454,7 @@ ntlGeometryObjModel::getTriangles(double t, vector *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 &verts, vector &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:"<getName()<<" frame:"<getName()<<" frame:"< *triangles, +void ParticleTracer::getTriangles(double time, vector *triangles, vector *vertices, vector *normals, int objectId ) { @@ -326,7 +326,7 @@ void ParticleTracer::getTriangles(double t, vector *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 *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,"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="<=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 "<6)) { + //if(valmpIso->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"<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;j1.) 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:"<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 } -- cgit v1.2.3