diff options
author | Nils Thuerey <nils@thuerey.de> | 2006-03-29 11:35:54 +0400 |
---|---|---|
committer | Nils Thuerey <nils@thuerey.de> | 2006-03-29 11:35:54 +0400 |
commit | 0a63b3c0ca032d7bb26a97164f034a49499eee68 (patch) | |
tree | c14cbc74d3b0e6136dc060fe4e8a1bf16804a90b /intern | |
parent | 0d2902b1fe50cd27f0046a4c0ecf140e6e1284a8 (diff) |
Several minor fixes:
- Added part of Austin's msvc8 fixes (vector::erase function
was "misused"), hopefully compiles better now.
- Ctrl-b now also bakes a selected fluidsim domain
similar to the softbodies.
- Added surface smoothing option for domains: default is
1, higher values result in a smoother surface (and probably
slightly higher comupation times), while 0 means the surface
is not modified at all.
- Added BLENDER_ELBEEMBOBJABORT environment variable in readBobj,
if >0 quits blender when a not yet existing fluidsim
frame should be loaded. Useful for rendering simulations
as far as possible from the command line.
- Surface normals pointer is now set to NULL in readfile.c
- Fixed win32 error string handling, now uses a function
to return the string from the solver.
- Fixed fluidsim particle halo scaling problem.
- Solver update
Diffstat (limited to 'intern')
22 files changed, 518 insertions, 373 deletions
diff --git a/intern/elbeem/extern/LBM_fluidsim.h b/intern/elbeem/extern/LBM_fluidsim.h index b879dd41396..c93d1292e47 100644 --- a/intern/elbeem/extern/LBM_fluidsim.h +++ b/intern/elbeem/extern/LBM_fluidsim.h @@ -68,16 +68,6 @@ int performElbeemSimulation(char *cfgfilename); void fluidsimGetAxisAlignedBB(struct Mesh *mesh, float obmat[][4], /*RET*/ float start[3], /*RET*/ float size[3], /*RET*/ struct Mesh **bbmesh ); -// implemented in intern/elbeem/utilities.cpp -/* set elbeem debug output level (0=off to 10=full on) */ -void elbeemSetDebugLevel(int level); -/* elbeem debug output function */ -void elbeemDebugOut(char *msg); - -/* estimate how much memory a given setup will require */ -double elbeemEstimateMemreq(int res, - float sx, float sy, float sz, - int refine, char *retstr); #endif diff --git a/intern/elbeem/extern/elbeem.h b/intern/elbeem/extern/elbeem.h index 83a9cb99afc..f0d154c18a4 100644 --- a/intern/elbeem/extern/elbeem.h +++ b/intern/elbeem/extern/elbeem.h @@ -73,6 +73,10 @@ typedef struct elbeemSimulationSettings { /* global transformation to apply to fluidsim mesh */ float surfaceTrafo[4*4]; + + /* development variables, testing for upcoming releases...*/ + float farFieldSize; + } elbeemSimulationSettings; @@ -122,8 +126,12 @@ extern "C" { // reset elbeemSimulationSettings struct with defaults void elbeemResetSettings(struct elbeemSimulationSettings*); -// start fluidsim init +// start fluidsim init (returns !=0 upon failure) int elbeemInit(struct elbeemSimulationSettings*); +// get failure message during simulation or init +// if an error occured (the string is copied into buffer, +// max. length = 256 chars ) +void elbeemGetErrorString(char *buffer); // reset elbeemMesh struct with zeroes void elbeemResetMesh(struct elbeemMesh*); @@ -135,17 +143,36 @@ int elbeemAddMesh(struct elbeemMesh*); int elbeemSimulate(void); -// helper function - simplify animation channels +// helper functions + +// simplify animation channels // returns if the channel and its size changed int elbeemSimplifyChannelFloat(float *channel, int *size); int elbeemSimplifyChannelVec3(float *channel, int *size); +// helper functions implemented in utilities.cpp + +/* set elbeem debug output level (0=off to 10=full on) */ +void elbeemSetDebugLevel(int level); +/* elbeem debug output function, prints if debug level >0 */ +void elbeemDebugOut(char *msg); + +/* estimate how much memory a given setup will require */ +double elbeemEstimateMemreq(int res, + float sx, float sy, float sz, + int refine, char *retstr); + + + #ifdef __cplusplus } #endif // __cplusplus + + /******************************************************************************/ -// internal defines, do not use for setting up simulation +// internal defines, do not use for initializing elbeemMesh +// structs, for these use OB_xxx defines above /*! fluid geometry init types */ #define FGI_FLAGSTART 16 diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp index 8e96879737f..daee5f831d6 100644 --- a/intern/elbeem/intern/elbeem.cpp +++ b/intern/elbeem/intern/elbeem.cpp @@ -67,6 +67,8 @@ void elbeemResetSettings(elbeemSimulationSettings *set) { set->generateVertexVectors = 0; set->surfaceSmoothing = 1.; + set->farFieldSize = 0.; + // init identity for(int i=0; i<16; i++) set->surfaceTrafo[i] = 0.0; for(int i=0; i<4; i++) set->surfaceTrafo[i*4+i] = 1.0; @@ -87,6 +89,13 @@ int elbeemInit(elbeemSimulationSettings *settings) { return 0; } +// error message access +extern "C" +void elbeemGetErrorString(char *buffer) { + if(!buffer) return; + strncpy(buffer,gElbeemErrorString,256); +} + // reset elbeemMesh struct with zeroes extern "C" void elbeemResetMesh(elbeemMesh *mesh) { diff --git a/intern/elbeem/intern/elbeem.h b/intern/elbeem/intern/elbeem.h index 83a9cb99afc..f0d154c18a4 100644 --- a/intern/elbeem/intern/elbeem.h +++ b/intern/elbeem/intern/elbeem.h @@ -73,6 +73,10 @@ typedef struct elbeemSimulationSettings { /* global transformation to apply to fluidsim mesh */ float surfaceTrafo[4*4]; + + /* development variables, testing for upcoming releases...*/ + float farFieldSize; + } elbeemSimulationSettings; @@ -122,8 +126,12 @@ extern "C" { // reset elbeemSimulationSettings struct with defaults void elbeemResetSettings(struct elbeemSimulationSettings*); -// start fluidsim init +// start fluidsim init (returns !=0 upon failure) int elbeemInit(struct elbeemSimulationSettings*); +// get failure message during simulation or init +// if an error occured (the string is copied into buffer, +// max. length = 256 chars ) +void elbeemGetErrorString(char *buffer); // reset elbeemMesh struct with zeroes void elbeemResetMesh(struct elbeemMesh*); @@ -135,17 +143,36 @@ int elbeemAddMesh(struct elbeemMesh*); int elbeemSimulate(void); -// helper function - simplify animation channels +// helper functions + +// simplify animation channels // returns if the channel and its size changed int elbeemSimplifyChannelFloat(float *channel, int *size); int elbeemSimplifyChannelVec3(float *channel, int *size); +// helper functions implemented in utilities.cpp + +/* set elbeem debug output level (0=off to 10=full on) */ +void elbeemSetDebugLevel(int level); +/* elbeem debug output function, prints if debug level >0 */ +void elbeemDebugOut(char *msg); + +/* estimate how much memory a given setup will require */ +double elbeemEstimateMemreq(int res, + float sx, float sy, float sz, + int refine, char *retstr); + + + #ifdef __cplusplus } #endif // __cplusplus + + /******************************************************************************/ -// internal defines, do not use for setting up simulation +// internal defines, do not use for initializing elbeemMesh +// structs, for these use OB_xxx defines above /*! fluid geometry init types */ #define FGI_FLAGSTART 16 diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp index ed0fab289d8..07bc6b7c855 100644 --- a/intern/elbeem/intern/isosurface.cpp +++ b/intern/elbeem/intern/isosurface.cpp @@ -37,7 +37,7 @@ IsoSurface::IsoSurface(double iso) : mInitDone(false), mSmoothSurface(0.0), mSmoothNormals(0.0), mAcrossEdge(), mAdjacentFaces(), - mCutoff(-1), // off by default + mCutoff(-1), mCutArray(NULL),// off by default mFlagCnt(1), mSCrad1(0.), mSCrad2(0.), mSCcenter(0.) { @@ -152,6 +152,7 @@ void IsoSurface::triangulate( void ) const int cubieOffsetZ[8] = { 0,0,0,0, 1,1,1,1 }; + const int coAdd=2; // let the cubes march pz = mStart[2]-gsz*0.5; for(int k=1;k<(mSizez-2);k++) { @@ -259,12 +260,17 @@ void IsoSurface::triangulate( void ) } - const int coAdd=2; - if(i<coAdd+mCutoff) continue; - if(j<coAdd+mCutoff) continue; - if((mCutoff>0) && (k<coAdd)) continue; - if(i>mSizex-2-coAdd-mCutoff) continue; - if(j>mSizey-2-coAdd-mCutoff) continue; + if( (i<coAdd+mCutoff) || + (j<coAdd+mCutoff) || + ((mCutoff>0) && (k<coAdd)) ||// bottom layer + (i>mSizex-2-coAdd-mCutoff) || + (j>mSizey-2-coAdd-mCutoff) ) { + if(mCutArray) { + if(k < mCutArray[j*this->mSizex+i]) continue; + } else { + continue; + } + } // Create the triangles... for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { diff --git a/intern/elbeem/intern/isosurface.h b/intern/elbeem/intern/isosurface.h index 90a4e11e95f..5e5115ef236 100644 --- a/intern/elbeem/intern/isosurface.h +++ b/intern/elbeem/intern/isosurface.h @@ -95,6 +95,8 @@ class IsoSurface : //! cutoff border area int mCutoff; + //! cutoff heigh values + int *mCutArray; //! trimesh vars vector<int> flags; @@ -160,6 +162,8 @@ class IsoSurface : } //! set cut off border inline void setCutoff(int set) { mCutoff = set; }; + //! set cut off border + inline void setCutArray(int *set) { mCutArray = set; }; //! OpenGL viz "interface" unsigned int getIsoVertexCount() { diff --git a/intern/elbeem/intern/ntl_blenderdumper.cpp b/intern/elbeem/intern/ntl_blenderdumper.cpp index fce5a085f59..ecd1967c877 100644 --- a/intern/elbeem/intern/ntl_blenderdumper.cpp +++ b/intern/elbeem/intern/ntl_blenderdumper.cpp @@ -91,7 +91,7 @@ int ntlBlenderDumper::renderScene( void ) vector<ntlTriangle> Triangles; vector<ntlVec3Gfx> Vertices; vector<ntlVec3Gfx> VertNormals; - errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) ); + //errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) ); /* init geometry array, first all standard objects */ int idCnt = 0; // give IDs to objects @@ -107,13 +107,13 @@ int ntlBlenderDumper::renderScene( void ) if(tid & GEOCLASSTID_SHADER) { ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter); //dynamic_cast<ntlGeometryShader*>(*iter); hideObjs.push_back( (*iter)->getName() ); - geoshad->notifyShaderOfDump(glob->getAniCount(),nrStr,glob->getOutFilename()); + geoshad->notifyShaderOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename()); for (vector<ntlGeometryObject*>::iterator siter = geoshad->getObjectsBegin(); siter != geoshad->getObjectsEnd(); siter++) { if(debugOut) debMsgStd("ntlBlenderDumper::BuildScene",DM_MSG,"added shader geometry "<<(*siter)->getName(), 8); - (*siter)->notifyOfDump(glob->getAniCount(),nrStr,glob->getOutFilename()); + (*siter)->notifyOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename(), this->mSimulationTime); bool doDump = false; bool isPreview = false; // only dump final&preview surface meshes diff --git a/intern/elbeem/intern/ntl_geometryobject.cpp b/intern/elbeem/intern/ntl_geometryobject.cpp index 7bbae5d87a1..5ab13a946d3 100644 --- a/intern/elbeem/intern/ntl_geometryobject.cpp +++ b/intern/elbeem/intern/ntl_geometryobject.cpp @@ -157,9 +157,9 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob) /*! notify object that dump is in progress (e.g. for particles) */ // default action - do nothing... -void ntlGeometryObject::notifyOfDump(int frameNr,char *frameNrStr,string outfilename) { +void ntlGeometryObject::notifyOfDump(int dumtp, int frameNr,char *frameNrStr,string outfilename, double simtime) { bool debugOut=false; - if(debugOut) debMsgStd("ntlGeometryObject::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<","<<frameNr<<" to "<<outfilename, 10); // DEBUG + if(debugOut) debMsgStd("ntlGeometryObject::notifyOfDump",DM_MSG," dt:"<<dumtp<<" obj:"<<this->getName()<<" frame:"<<frameNrStr<<","<<frameNr<<",t"<<simtime<<" to "<<outfilename, 10); // DEBUG } /*****************************************************************************/ diff --git a/intern/elbeem/intern/ntl_geometryobject.h b/intern/elbeem/intern/ntl_geometryobject.h index 157d160c0e4..ae81d748818 100644 --- a/intern/elbeem/intern/ntl_geometryobject.h +++ b/intern/elbeem/intern/ntl_geometryobject.h @@ -16,6 +16,8 @@ class ntlRenderGlobals; class ntlTriangle; +#define DUMP_FULLGEOMETRY 1 +#define DUMP_PARTIAL 2 class ntlGeometryObject : public ntlGeometryClass { @@ -38,7 +40,7 @@ class ntlGeometryObject : public ntlGeometryClass vector<ntlVec3Gfx> *normals, int objectId ) = 0; /*! notify object that dump is in progress (e.g. for particles) */ - virtual void notifyOfDump(int frameNr,char *frameNrStr,string outfilename); + virtual void notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime); /*! Search the material for this object from the material list */ void searchMaterial(vector<ntlMaterial *> *mat); diff --git a/intern/elbeem/intern/ntl_geometryshader.h b/intern/elbeem/intern/ntl_geometryshader.h index 7d767b0f8b9..6adc5629a2d 100644 --- a/intern/elbeem/intern/ntl_geometryshader.h +++ b/intern/elbeem/intern/ntl_geometryshader.h @@ -40,7 +40,7 @@ class ntlGeometryShader : virtual vector<ntlGeometryObject *>::iterator getObjectsEnd() { return mObjects.end(); } /*! notify object that dump is in progress (e.g. for field dump) */ - virtual void notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename) = 0; + virtual void notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0; protected: diff --git a/intern/elbeem/intern/ntl_world.h b/intern/elbeem/intern/ntl_world.h index 2256e85ec96..131f1828601 100644 --- a/intern/elbeem/intern/ntl_world.h +++ b/intern/elbeem/intern/ntl_world.h @@ -380,9 +380,6 @@ private: bool mSingleFrameMode; //! filename for single frame mode string mSingleFrameFilename; - - /*! Two random number streams for photon generation (one for the directions, the other for russion roulette) */ - //ntlRandomStream *mpRndDirections, *mpRndRoulette; }; diff --git a/intern/elbeem/intern/particletracer.cpp b/intern/elbeem/intern/particletracer.cpp index 0a5af90564a..828674675e2 100644 --- a/intern/elbeem/intern/particletracer.cpp +++ b/intern/elbeem/intern/particletracer.cpp @@ -20,6 +20,8 @@ #include <zlib.h> +// particle object id counter +int ParticleObjectIdCnt = 1; /****************************************************************************** * Standard constructor @@ -34,19 +36,21 @@ ParticleTracer::ParticleTracer() : mPartScale(1.0) , mPartHeadDist( 0.5 ), mPartTailDist( -4.5 ), mPartSegments( 4 ), mValueScale(0), mValueCutoffTop(0.0), mValueCutoffBottom(0.0), - mDumpParts(0), mShowOnly(0), mpTrafo(NULL) + mDumpParts(0), mDumpText(0), mDumpTextFile(""), mShowOnly(0), + mNumInitialParts(0), mpTrafo(NULL) { // check env var -#ifdef ELBEEM_PLUGIN - mDumpParts=1; // default on -#endif // ELBEEM_PLUGIN - if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG! - int set = atoi( getenv("ELBEEM_DUMPPARTICLE") ); - if((set>=0)&&(set!=mDumpParts)) { - mDumpParts=set; - debMsgStd("ParticleTracer::notifyOfDump",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8); - } - } +//#ifdef ELBEEM_PLUGIN + //mDumpParts=1; // default on +//#endif // ELBEEM_PLUGIN + //if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG! + //int set = atoi( getenv("ELBEEM_DUMPPARTICLE") ); + //if((set>=0)&&(set!=mDumpParts)) { + //mDumpParts=set; + //debMsgStd("ParticleTrace",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8); + //} + //} + debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10); }; ParticleTracer::~ParticleTracer() { @@ -60,13 +64,9 @@ void ParticleTracer::parseAttrList(AttributeList *att) { AttributeList *tempAtt = mpAttrs; mpAttrs = att; - int mNumParticles =0; // UNUSED - int mTrailLength = 0; // UNUSED - int mTrailInterval= 0; // UNUSED - mNumParticles = mpAttrs->readInt("particles",mNumParticles, "ParticleTracer","mNumParticles", false); - mTrailLength = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false); - mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false); + mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false); + errMsg(" NUMP"," "<<mNumInitialParts); mPartScale = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false); mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false); mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false); @@ -75,18 +75,23 @@ void ParticleTracer::parseAttrList(AttributeList *att) mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false); mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false); - mDumpParts = mpAttrs->readInt ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false); + mDumpParts = mpAttrs->readInt ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false); + mDumpText = mpAttrs->readInt ("part_textdump",mDumpText, "ParticleTracer","mDumpText", false); mShowOnly = mpAttrs->readInt ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false); + mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false); string matPart; matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false); setMaterialName( matPart ); - // trail length has to be at least one, if anything should be displayed - //if((mNumParticles>0)&&(mTrailLength<2)) mTrailLength = 2; + + // unused... + int mTrailLength = 0; // UNUSED + int mTrailInterval= 0; // UNUSED + mTrailLength = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false); + mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false); // restore old list mpAttrs = tempAtt; - //mParts.resize(mTrailLength*mTrailInterval); } /****************************************************************************** @@ -146,7 +151,7 @@ void ParticleTracer::addParticle(float x, float y, float z) void ParticleTracer::cleanup() { // cleanup int last = (int)mParts.size()-1; - //for(vector<ParticleObject>::iterator pit= getParticlesBegin();pit!= getParticlesEnd(); pit++) { + if(mDumpText>0) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; } for(int i=0; i<=last; i++) { if( mParts[i].getActive()==false ) { @@ -158,47 +163,18 @@ void ParticleTracer::cleanup() { } /****************************************************************************** - * save particle positions before adding a new timestep - * copy "one index up", newest has to remain unmodified, it will be - * advanced after the next smiulation step + *! dump particles if desired *****************************************************************************/ -void ParticleTracer::savePreviousPositions() -{ - //debugOut(" PARTS SIZE "<<mParts.size() ,10); - /* - if(mTrailIntervalCounter==0) { - //errMsg("spp"," PARTS SIZE "<<mParts.size() ); - for(size_t l=mParts.size()-1; l>0; l--) { - if( mParts[l].size() != mParts[l-1].size() ) { - errFatal("ParticleTracer::savePreviousPositions","Invalid array sizes ["<<l<<"]="<<mParts[l].size()<< - " ["<<(l+1)<<"]="<<mParts[l+1].size() <<" , total "<< mParts.size() , SIMWORLD_GENERICERROR); - return; - } - - for(size_t i=0; i<mParts[l].size(); i++) { - mParts[l][i] = mParts[l-1][i]; - } - - } - } - mTrailIntervalCounter++; - if(mTrailIntervalCounter>=mTrailInterval) mTrailIntervalCounter = 0; - UNUSED!? */ -} - - -/*! dump particles if desired */ -void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilename) { +void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) { debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr, 10); // DEBUG - if(mDumpParts>0) { + if( + (dumptype==DUMP_FULLGEOMETRY)&& + (mDumpParts>0)) { // dump to binary file std::ostringstream boutfilename(""); - //boutfilename << ecrpath.str() << glob->getOutFilename() <<"_"<< this->getName() <<"_" << frameNrStr << ".obj"; - //boutfilename << outfilename <<"_particles_"<< this->getName() <<"_" << frameNrStr<< ".gz"; boutfilename << outfilename <<"_particles_" << frameNrStr<< ".gz"; - debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() - <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7); + debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7); // output to zipped file gzFile gzf; @@ -210,21 +186,26 @@ void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilenam numParts = 0; for(size_t i=0; i<mParts.size(); i++) { if(!mParts[i].getActive()) continue; + //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK numParts++; } gzwrite(gzf, &numParts, sizeof(numParts)); for(size_t i=0; i<mParts.size(); i++) { if(!mParts[i].getActive()) { continue; } - if(mParts[i].getLifeTime()<30) { continue; } //? CHECK + //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK ParticleObject *p = &mParts[i]; - int type = p->getType(); + //int type = p->getType(); // export whole type info + int type = p->getFlags(); // debug export whole type & status info ntlVec3Gfx pos = p->getPos(); float size = p->getSize(); if(type&PART_FLOAT) { // WARNING same handling for dump! // add one gridcell offset //pos[2] += 1.0; - } + } + // display as drop for now externally + //else if(type&PART_TRACER) { type |= PART_DROP; } + pos = (*mpTrafo) * pos; ntlVec3Gfx v = p->getVel(); @@ -240,6 +221,64 @@ void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilenam gzclose( gzf ); } } // dump? + + // dfor partial & full dump + if(mDumpText>0) { + // dump to binary file + std::ostringstream boutfilename(""); + if(mDumpTextFile.length()>1) { + boutfilename << mDumpTextFile << ".cpart2"; + } else { + boutfilename << outfilename <<"_particles" << ".cpart2"; + } + debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7); + + int numParts = 0; + // only dump bubble particles + for(size_t i=0; i<mParts.size(); i++) { + //if(!mParts[i].getActive()) continue; + //if(!(mParts[i].getType()&PART_BUBBLE)) continue; + numParts++; + } + + // output to text file + gzFile gzf; + if(frameNr==0) { + gzf = gzopen(boutfilename.str().c_str(), "w0"); + + gzprintf( gzf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts); + // fixed time scale for now + gzprintf( gzf, "T %f \n\n", 1.0); + } else { + gzf = gzopen(boutfilename.str().c_str(), "a+0"); + } + + // add current set + if(gzf) { + gzprintf( gzf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", frameNr, simtime, numParts ); + gzprintf( gzf, "S %f \n\n", simtime ); + + for(size_t i=0; i<mParts.size(); i++) { + ParticleObject *p = &mParts[i]; + ntlVec3Gfx pos = p->getPos(); + float size = p->getSize(); + if(!mParts[i].getActive()) { size=0.; } // switch "off" + + pos = (*mpTrafo) * pos; + ntlVec3Gfx v = p->getVel(); + v[0] *= mpTrafo->value[0][0]; + v[1] *= mpTrafo->value[1][1]; + v[2] *= mpTrafo->value[2][2]; + + gzprintf( gzf, "P %f %f %f \n", pos[0],pos[1],pos[2] ); + gzprintf( gzf, "s %f \n", size ); + gzprintf( gzf, "\n", size ); + } + + gzprintf( gzf, "# %d end ", frameNr ); + gzclose( gzf ); + } + } } @@ -293,6 +332,7 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles, case 2: if(!(type&PART_DROP)) continue; break; case 3: if(!(type&PART_INTER)) continue; break; case 4: if(!(type&PART_FLOAT)) continue; break; + case 5: if(!(type&PART_TRACER)) continue; break; } } else { // by default dont display inter diff --git a/intern/elbeem/intern/particletracer.h b/intern/elbeem/intern/particletracer.h index b9179b68dc0..453c1a793c0 100644 --- a/intern/elbeem/intern/particletracer.h +++ b/intern/elbeem/intern/particletracer.h @@ -16,24 +16,31 @@ template<class Scalar> class ntlMatrix4x4; #define PART_DROP (1<< 2) #define PART_INTER (1<< 3) #define PART_FLOAT (1<< 4) +#define PART_TRACER (1<< 5) // particle state #define PART_IN (1<< 8) #define PART_OUT (1<< 9) #define PART_INACTIVE (1<<10) +// defines for particle movement +#define MOVE_FLOATS 1 +#define FLOAT_JITTER 0.03 + +extern int ParticleObjectIdCnt; + //! A single particle class ParticleObject { public: //! Standard constructor inline ParticleObject(ntlVec3Gfx mp) : - mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0) { }; + mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0) { mId = ParticleObjectIdCnt++; }; //! Copy constructor inline ParticleObject(const ParticleObject &a) : mPos(a.mPos), mVel(a.mVel), mSize(a.mSize), mStatus(a.mStatus), - mLifeTime(a.mLifeTime) { }; + mLifeTime(a.mLifeTime) { mId = ParticleObjectIdCnt++; }; //! Destructor inline ~ParticleObject() { /* empty */ }; @@ -81,8 +88,12 @@ class ParticleObject //! set type (lower byte) inline void setLifeTime(int set) { mLifeTime = set; } + inline int getId() const { return mId; } + protected: + /*! only for debugging */ + int mId; /*! the particle position */ ntlVec3Gfx mPos; /*! the particle velocity */ @@ -109,9 +120,6 @@ class ParticleTracer : //! add a particle at this position void addParticle(float x, float y, float z); - //! save particle positions before adding a new timestep - void savePreviousPositions(); - //! draw the particle array void draw(); @@ -125,6 +133,9 @@ class ParticleTracer : //! get the number of particles inline int getNumParticles() { return mParts.size(); } + //! set/get the number of particles + inline void setNumInitialParticles(int set) { mNumInitialParts=set; } + inline int getNumInitialParticles() { return mNumInitialParts; } //! iterate over all newest particles (for advancing positions) inline vector<ParticleObject>::iterator getParticlesBegin() { return mParts.begin(); } @@ -149,7 +160,13 @@ class ParticleTracer : /*! set/get dump flag */ inline void setDumpParts(bool set) { mDumpParts = set; } - inline bool getDumpParts() { return mDumpParts; } + inline bool getDumpParts() { return mDumpParts; } + /*! set/get dump flag */ + inline void setDumpText(bool set) { mDumpText = set; } + inline bool getDumpText() { return mDumpText; } + /*! set/get dump text file */ + inline void setDumpTextFile(std::string set) { mDumpTextFile = set; } + inline std::string getDumpTextFile() { return mDumpTextFile; } //! set the particle scaling factor inline void setPartScale(float set) { mPartScale = set; } @@ -161,7 +178,7 @@ class ParticleTracer : vector<ntlVec3Gfx> *vertices, vector<ntlVec3Gfx> *normals, int objectId ); - virtual void notifyOfDump(int frameNr,char *frameNrStr,string outfilename); + virtual void notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename,double simtime); // free deleted particles void cleanup(); @@ -194,8 +211,14 @@ class ParticleTracer : /*! dump particles (or certain types of) to disk? */ int mDumpParts; + /*! dump particles (or certain types of) to disk? */ + int mDumpText; + /*! text dump output file */ + std::string mDumpTextFile; /*! show only a certain type (debugging) */ int mShowOnly; + /*! no. of particles to init */ + int mNumInitialParts; //! transform matrix ntlMatrix4x4<gfxReal> *mpTrafo; diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp index 481eb836ad6..116d68d6824 100644 --- a/intern/elbeem/intern/simulation_object.cpp +++ b/intern/elbeem/intern/simulation_object.cpp @@ -156,6 +156,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) mpLbm->setGeoEnd( mGeoEnd ); mpLbm->setRenderGlobals( mpGlob ); mpLbm->setName( getName() + "_lbm" ); + mpLbm->setParticleTracer( mpParts ); if(mpElbeemSettings) { // set further settings from API struct init mpLbm->setSmoothing(1.0 * mpElbeemSettings->surfaceSmoothing, 1.0 * mpElbeemSettings->surfaceSmoothing); @@ -173,6 +174,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) mpLbm->setDomainBound(dinitType); mpLbm->setDomainPartSlip(mpElbeemSettings->obstaclePartslip); mpLbm->setDumpVelocities(mpElbeemSettings->generateVertexVectors); + mpLbm->setFarFieldSize(mpElbeemSettings->farFieldSize); debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<<dinitType<<" ps="<<mpElbeemSettings->obstaclePartslip<<" vv"<<mpElbeemSettings->generateVertexVectors<<","<<mpLbm->getDumpVelocities(), 9 ); debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<<mpLbm->getGenerateParticles(),10); @@ -236,7 +238,6 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) mpParts->setCastShadows( false ); mpParts->setReceiveShadows( false ); mpParts->searchMaterial( glob->getMaterials() ); - mpLbm->initParticles(mpParts); // this has to be inited here - before, the values might be unknown ntlGeometryObject *surf = mpLbm->getSurfaceGeoObj(); @@ -288,9 +289,6 @@ void SimulationObject::step( void ) if(mpParam->getCurrentAniFrameTime()>0.0) { // dont advance for stopped time mpLbm->step(); - - mpParts->savePreviousPositions(); - mpLbm->advanceParticles(mpParts); mTime += mpParam->getTimestep(); } if(mpLbm->getPanic()) mPanic = true; @@ -399,8 +397,8 @@ void SimulationObject::setMouseClick() } /*! notify object that dump is in progress (e.g. for field dump) */ -void SimulationObject::notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename) { +void SimulationObject::notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) { if(!mpLbm) return; - mpLbm->notifySolverOfDump(frameNr,frameNrStr,outfilename); + mpLbm->notifySolverOfDump(dumptype, frameNr,frameNrStr,outfilename); } diff --git a/intern/elbeem/intern/simulation_object.h b/intern/elbeem/intern/simulation_object.h index 4e7bc19294d..b8d1e90596d 100644 --- a/intern/elbeem/intern/simulation_object.h +++ b/intern/elbeem/intern/simulation_object.h @@ -91,7 +91,7 @@ class SimulationObject : virtual int postGeoConstrInit(ntlRenderGlobals *glob) { return initializeLbmSimulation(glob); }; virtual int initializeShader() { /* ... */ return true; }; /*! notify object that dump is in progress (e.g. for field dump) */ - virtual void notifyShaderOfDump(int frameNr,char *frameNrStr,string outfilename); + virtual void notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename); /*! simluation interface: draw the simulation with OpenGL */ virtual void draw( void ) {}; virtual vector<ntlGeometryObject *>::iterator getObjectsBegin(); diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h index 7217964e7ba..b24183411ef 100644 --- a/intern/elbeem/intern/solver_class.h +++ b/intern/elbeem/intern/solver_class.h @@ -211,7 +211,7 @@ class LbmFsgrSolver : //! finish the init with config file values (allocate arrays...) virtual bool initializeSolver(); //! notify object that dump is in progress (e.g. for field dump) - virtual void notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename); + virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename); #if LBM_USE_GUI==1 //! show simulation info (implement LbmSolverInterface pure virtual func) @@ -267,9 +267,9 @@ class LbmFsgrSolver : /* simulation object interface, just calls stepMain */ virtual void step(); /*! init particle positions */ - virtual int initParticles(ParticleTracer *partt); + virtual int initParticles(); /*! move all particles */ - virtual void advanceParticles(ParticleTracer *partt ); + virtual void advanceParticles(); /*! debug object display (used e.g. for preview surface) */ @@ -436,8 +436,6 @@ class LbmFsgrSolver : int mForceTadapRefine; //! border cutoff value int mCutoff; - //! store particle tracer - ParticleTracer *mpParticles; // strict debug interface # if FSGR_STRICT_DEBUG==1 @@ -460,11 +458,10 @@ class LbmFsgrSolver : void initTestdata(); void destroyTestdata(); void handleTestdata(); - void exportTestdata(); void set3dHeight(int ,int ); public: // needed from testdata - void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat &retux, LbmFloat &retuy); + void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy); #endif // LBM_INCLUDE_TESTSOLVERS==1 public: // former LbmModelLBGK functions @@ -700,9 +697,8 @@ class LbmFsgrSolver : #define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set) ] #define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set) ] -// array data layouts -// standard array layout ----------------------------------------------------------------------------------------------- -#define ALSTRING "Standard Array Layout" +// array handling ----------------------------------------------------------------------------------------------- + //#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) ) #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) ) #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)]) diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp index abe9ef6994a..a417c33caa9 100644 --- a/intern/elbeem/intern/solver_init.cpp +++ b/intern/elbeem/intern/solver_init.cpp @@ -396,8 +396,6 @@ LbmFsgrSolver::LbmFsgrSolver() : mGaussw[n] = mGaussw[n]/totGaussw; } - mpParticles = NULL; - //addDrop(false,0,0); } /*****************************************************************************/ @@ -486,6 +484,8 @@ void LbmFsgrSolver::parseAttrList() #else // LBM_INCLUDE_TESTSOLVERS!=1 // off by default mUseTestdata = 0; + if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check + if(mUseTestdata) { mMaxRefine=0; } // force fsgr off #endif // LBM_INCLUDE_TESTSOLVERS!=1 } @@ -573,7 +573,14 @@ void LbmFsgrSolver::initLevelOmegas() *****************************************************************************/ bool LbmFsgrSolver::initializeSolver() { -// debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... (Layout:"<<ALSTRING<<") "<<this->mInitDone<<" "<<((int)this),1); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<this->mInitDone<<" "<<(void*)this,1); + + // init cppf stage + if(mCppfStage>0) { + this->mSizex *= mCppfStage; + this->mSizey *= mCppfStage; + this->mSizez *= mCppfStage; + } // size inits to force cubic cells and mult4 level dimensions // and make sure we dont allocate too much... @@ -584,9 +591,6 @@ bool LbmFsgrSolver::initializeSolver() double sizeReduction = 1.0; double memEstFromFunc = -1.0; string memreqStr(""); -#if LBM_INCLUDE_TESTSOLVERS==1 - if(mUseTestdata) { mMaxRefine=0; } // force fsgr off -#endif while(!memOk) { initGridSizes( this->mSizex, this->mSizey, this->mSizez, this->mvGeoStart, this->mvGeoEnd, mMaxRefine, PARALLEL); @@ -750,16 +754,25 @@ bool LbmFsgrSolver::initializeSolver() mMaxTimestep = this->mpParam->getTimestep(); // init isosurf + this->mpIso->setIsolevel( this->mIsoValue ); #if LBM_INCLUDE_TESTSOLVERS==1 if(mUseTestdata) { mpTest->setMaterialName( this->mpIso->getMaterialName() ); delete this->mpIso; this->mpIso = mpTest; + if(mpTest->mDebugvalue1>0.0) { // 3d off + mpTest->setIsolevel(-100.0); + } else { + mpTest->setIsolevel( this->mIsoValue ); + } } #endif // ELBEEM_PLUGIN!=1 - this->mpIso->setIsolevel( this->mIsoValue ); // approximate feature size with mesh resolution float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5; + // smooth vars defined in solver_interface, set by simulation object + // reset for invalid values... + if((this->mSmoothSurface<0.)||(this->mSmoothSurface>50.)) this->mSmoothSurface = 1.; + if((this->mSmoothNormals<0.)||(this->mSmoothNormals>50.)) this->mSmoothNormals = 1.; this->mpIso->setSmoothSurface( this->mSmoothSurface * featureSize ); this->mpIso->setSmoothNormals( this->mSmoothNormals * featureSize ); @@ -984,6 +997,12 @@ bool LbmFsgrSolver::initializeSolver() errMsg("LbmFsgrSolver::init","No preview in 2D allowed!"); this->mOutputSurfacePreview = 0; } } +#if LBM_USE_GUI==1 + if(this->mOutputSurfacePreview) { + errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0"); + this->mOutputSurfacePreview = 0; } +#endif // LBM_USE_GUI==1 + if(this->mOutputSurfacePreview) { // same as normal one, but use reduced size @@ -1017,13 +1036,17 @@ bool LbmFsgrSolver::initializeSolver() // now really done... - debugOut("LbmFsgrSolver::initialize : Init done ...",10); + debugOut("LbmFsgrSolver::initialize : SurfaceGen: SmsOrg("<<this->mSmoothSurface<<","<<this->mSmoothNormals<<","<<featureSize<<"), Iso("<<this->mpIso->getSmoothSurface()<<","<<this->mpIso->getSmoothNormals()<<") ",10); + debugOut("LbmFsgrSolver::initialize : Init done ... ",10); this->mInitDone = 1; #if LBM_INCLUDE_TESTSOLVERS==1 initTestdata(); #endif // ELBEEM_PLUGIN!=1 + // not inited? dont use... + if(mCutoff<0) mCutoff=0; + initParticles(); return true; } diff --git a/intern/elbeem/intern/solver_interface.cpp b/intern/elbeem/intern/solver_interface.cpp index 04e0456286b..ee2f8f281be 100644 --- a/intern/elbeem/intern/solver_interface.cpp +++ b/intern/elbeem/intern/solver_interface.cpp @@ -34,7 +34,7 @@ LbmSolverInterface::LbmSolverInterface() : mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop( (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ), mInitDone( false ), mInitDensityGradient( false ), - mpAttrs( NULL ), mpParam( NULL ), + mpAttrs( NULL ), mpParam( NULL ), mpParticles(NULL), mNumParticlesLost(0), mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0), mDebugVelScale( 0.01 ), mNodeInfoString("+"), @@ -53,7 +53,7 @@ LbmSolverInterface::LbmSolverInterface() : mDumpVelocities(false), mMarkedCells(), mMarkedCellIndex(0), mDomainBound("noslip"), mDomainPartSlipValue(0.1), - mTForceStrength(0.0) + mTForceStrength(0.0), mFarFieldSize(0.), mCppfStage(0) { #if ELBEEM_PLUGIN==1 if(gDebugLevel<=1) mSilent = true; @@ -242,7 +242,12 @@ void LbmSolverInterface::parseStdAttrList() { // new test vars mTForceStrength = mpAttrs->readFloat("tforcestrength", mTForceStrength,"LbmSolverInterface", "mTForceStrength", false); + mFarFieldSize = mpAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false); + // old compat + float sizeScale = mpAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false); + if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); } + mCppfStage = mpAttrs->readFloat("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false); mPartGenProb = mpAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false); } diff --git a/intern/elbeem/intern/solver_interface.h b/intern/elbeem/intern/solver_interface.h index c86bbce2a89..ee4ecd79077 100644 --- a/intern/elbeem/intern/solver_interface.h +++ b/intern/elbeem/intern/solver_interface.h @@ -260,7 +260,7 @@ class LbmSolverInterface virtual bool initializeSolver() =0; /*! notify object that dump is in progress (e.g. for field dump) */ - virtual void notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename) = 0; + virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0; /*! parse a boundary flag string */ CellFlagType readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed); @@ -273,8 +273,8 @@ class LbmSolverInterface virtual void prepareVisualization() { /* by default off */ }; /*! particle handling */ - virtual int initParticles(ParticleTracer *partt) = 0; - virtual void advanceParticles(ParticleTracer *partt ) = 0; + virtual int initParticles() = 0; + virtual void advanceParticles() = 0; /*! get surface object (NULL if no surface) */ ntlGeometryObject* getSurfaceGeoObj() { return mpIso; } @@ -329,6 +329,9 @@ class LbmSolverInterface inline void setParametrizer(Parametrizer *set) { mpParam = set; } /*! get parametrizer pointer */ inline Parametrizer *getParametrizer() { return mpParam; } + /*! get/set particle pointer */ + inline void setParticleTracer(ParticleTracer *set) { mpParticles = set; } + inline ParticleTracer *getParticleTracer() { return mpParticles; } /*! set density gradient init from e.g. init test cases */ inline void setInitDensityGradient(bool set) { mInitDensityGradient = set; } @@ -379,6 +382,12 @@ class LbmSolverInterface //! set/get dump velocities flag inline void setDomainPartSlip(LbmFloat set) { mDomainPartSlipValue = set; } inline LbmFloat getDomainPartSlip() const { return mDomainPartSlipValue; } + //! set/get far field size + inline void setFarFieldSize(LbmFloat set) { mFarFieldSize = set; } + inline LbmFloat getFarFieldSize() const { return mFarFieldSize; } + //! set/get cp stage + inline void setCpStage(int set) { mCppfStage = set; } + inline int getCpStage() const { return mCppfStage; } // cell iterator interface @@ -474,6 +483,8 @@ class LbmSolverInterface /*! get parameters from this parametrize in finishInit */ Parametrizer *mpParam; + //! store particle tracer + ParticleTracer *mpParticles; /*! number of particles lost so far */ int mNumParticlesLost; @@ -553,6 +564,10 @@ class LbmSolverInterface //! test vars // strength of applied force LbmFloat mTForceStrength; + // + LbmFloat mFarFieldSize; + // + int mCppfStage; }; diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp index c8f0678a6f0..10fca60c546 100644 --- a/intern/elbeem/intern/solver_main.cpp +++ b/intern/elbeem/intern/solver_main.cpp @@ -108,7 +108,7 @@ void LbmFsgrSolver::stepMain() int avgcls = (int)(mAvgNumUsedCells/(LONGINT)this->mStepCnt); debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime<< //debMsgDirect( - "mlsups(curr:"<<this->mMLSUPS<< + " mlsups(curr:"<<this->mMLSUPS<< " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<< " totcls:"<<(this->mNumUsedCells)<< sepStr<< " avgcls:"<< avgcls<< sepStr<< @@ -122,9 +122,6 @@ void LbmFsgrSolver::stepMain() " probs:"<<mNumProblems<< sepStr<< " simt:"<<mSimulationTime<< sepStr<< " for '"<<this->mName<<"' " , 10); - - debMsgDirect(std::endl); - debMsgDirect(this->mStepCnt<<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<") "); debMsgDirect(std::endl); // nicer output @@ -208,6 +205,9 @@ void LbmFsgrSolver::stepMain() if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); } #endif + // advance positions with current grid + advanceParticles(); + // one of the last things to do - adapt timestep // was in fineAdvance before... if(mTimeAdap) { @@ -220,6 +220,13 @@ void LbmFsgrSolver::stepMain() if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; } if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; } #endif // WIN32 + + // output total step time + timeend = getTime(); + debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt + <<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), " + <<" totst:"<<getTimeString(timeend-timestart), 3); + //#endif // ELBEEM_PLUGIN!=1 } @@ -283,11 +290,12 @@ LbmFsgrSolver::mainLoop(int lev) int calcCellsFilled = this->mNumFilledCells; int calcCellsEmptied = this->mNumEmptiedCells; int calcNumUsedCells = this->mNumUsedCells; + const int cutMin = 1; + const int cutConst = mCutoff+1; # if LBM_INCLUDE_TESTSOLVERS==1 - if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { - // 3d region off... quit - this->mpIso->setIsolevel(-100.0); return; } + // 3d region off... quit + if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; } #endif // ELBEEM_PLUGIN!=1 //printLbmCell(lev, 6,6,16, mLevel[lev].setCurr ); // DEBUG @@ -814,13 +822,18 @@ LbmFsgrSolver::mainLoop(int lev) // for inflow no pargen test // NOBUBBB! if((this->mInitDone) //&&(mUseTestdata) - && (!((oldFlag|newFlag)&CFNoNbEmpty)) + //&& (!((oldFlag|newFlag)&CFNoNbEmpty)) && (!((oldFlag|newFlag)&CFNoDelete)) && (this->mPartGenProb>0.0)) { + bool doAdd = true; bool bndOk=true; - if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)|| - (j<1+mCutoff)||(j>this->mSizey-1-mCutoff)|| - (k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; } + //if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)|| + //(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)|| + //(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; } + if( (i<cutMin)||(i>this->mSizex-cutMin)|| + (j<cutMin)||(j>this->mSizey-cutMin)|| + (k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; } + if(!bndOk) doAdd=false; LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep); LbmFloat rux = (ux * realWorldFac); @@ -828,12 +841,19 @@ LbmFsgrSolver::mainLoop(int lev) LbmFloat ruz = (uz * realWorldFac); LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz)); // WHMOD - const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor! - bool doAdd = true; LbmFloat prob = (rand()/(RAND_MAX+1.0)); LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl; - if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>2.5) ) { + + // reduce probability in outer region + if( (i<cutConst)||(i>this->mSizex-cutConst) ){ prob *= 0.5; } + if( (j<cutConst)||(j>this->mSizey-cutConst) ){ prob *= 0.5; } + if( (k<cutConst)||(k>this->mSizez-cutConst) ){ prob *= 0.5; } + +//#define RWVEL_THRESH 1.0 +#define RWVEL_THRESH 1.5 +#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5) + if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { // add } else { doAdd = false; // dont... @@ -851,12 +871,12 @@ LbmFsgrSolver::mainLoop(int lev) // "wind" disturbance // use realworld relative velocity here instead? - if( - ((rl>1.0) && (lcsmqo<P_LCSMQO)) // normal checks - ||(k>this->mSizez-SLOWDOWNREGION) ) { + if( (doAdd) && ( + ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) // normal checks + ||(k>this->mSizez-SLOWDOWNREGION) )) { LbmFloat nuz = uz; LbmFloat jdf; // = 0.05 * (rand()/(RAND_MAX+1.0)); - if(rl>1.0) jdf *= (rl-1.0); + if(rl>RWVEL_WINDTHRESH) jdf *= (rl-RWVEL_WINDTHRESH); if(k>this->mSizez-SLOWDOWNREGION) { // special case LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) ); @@ -881,7 +901,6 @@ LbmFsgrSolver::mainLoop(int lev) if(usqr<0.0001) doAdd=false; // TODO check!? // if outside, and 20% above sea level, delete, TODO really check level? //if((!bndOk)&&((LbmFloat)k>pTest->mFluidHeight*1.5)) { doAdd=true; } // FORCEDISSOLVE - //? if(!bndOk) doAdd=false; //if(this->mStepCnt>700) errMsg("DFJITT"," at "<<PRINT_IJK<<"rwl:"<<rl<<" usqr:"<<usqr <<" qo:"<<lcsmqo<<" add="<<doAdd ); if( (doAdd) ) { // ADD DROP @@ -911,15 +930,14 @@ LbmFsgrSolver::mainLoop(int lev) LbmFloat srcj = j+0.5+jpy; // TEST? + (pv[1]*1.41); LbmFloat srck = k+0.5+jpz; // TEST? + (pv[2]*1.41); int type=0; - //if((s%3)!=2) { - - type=PART_DROP; - // drop - srci += (pv[0]*1.41); - srcj += (pv[1]*1.41); - srck += (pv[2]*1.41); - if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction - //} else { type=PART_FLOAT; } + //if((s%3)!=2) {} else { type=PART_FLOAT; } + //type = PART_DROP; + type = PART_INTER; + // drop + /*srci += (pv[0]*1.41); + srcj += (pv[1]*1.41); + srck += (pv[2]*1.41); + if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction */ pv *= len; LbmFloat size = 1.0+ 9.0* (rand()/(RAND_MAX+1.0)); @@ -929,18 +947,20 @@ LbmFsgrSolver::mainLoop(int lev) //? mpParticles->getLast()->advanceVel(); // advance a bit outwards mpParticles->getLast()->setStatus(PART_IN); mpParticles->getLast()->setType(type); - //mpParticles->getLast()->setType(PART_INTER); //if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT); mpParticles->getLast()->setSize(size); //errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<PRINT_VEC(ux,uy,uz) <<" RWu="<<PRINT_VEC(rux,ruy,ruz)<<" add"<<doAdd<<" pvel="<<pv ); //errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) ); //if(!bndOk){ mass -= val2fac*size*0.02; } // FORCEDISSOLVE - mass -= val2fac*size*0.0015; // NTEST! + //const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor! + //mass -= val2fac*size*0.0015; // NTEST! //mass -= val2fac*size*0.001; // NTEST! + mass -= size*0.0020; // NTEST! #if LBMDIM==2 mpParticles->getLast()->setVel(pv[0],pv[1],0.0); mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5)); #endif // LBMDIM==2 + //errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel() ); } // multiple parts } // doAdd } // */ @@ -2487,17 +2507,19 @@ void LbmFsgrSolver::reinitFlags( int workSet ) { /* remove empty interface cells that are not allowed to be removed anyway * this is important, otherwise the dreaded cell-type-flickering can occur! */ - for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); - iter != mListEmpty.end(); iter++ ) { - int i=iter->x, j=iter->y, k=iter->z; + //for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); iter != mListEmpty.end(); iter++ ) { + //int i=iter->x, j=iter->y, k=iter->z; + //iter = mListEmpty.erase(iter); iter--; // and continue with next... + for(int n=0; n<(int)mListEmpty.size(); n++) { + int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z; if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) { - // remove entry - if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM - iter = mListEmpty.erase(iter); - iter--; // and continue with next... - // treat as "new inter" addToNewInterList(i,j,k); + // remove entry + if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM + if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; + mListEmpty.pop_back(); + n--; } } diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h index 589f2197422..605bc02c0cf 100644 --- a/intern/elbeem/intern/solver_relax.h +++ b/intern/elbeem/intern/solver_relax.h @@ -17,102 +17,9 @@ -#if LBM_INCLUDE_TESTSOLVERS!=1 - // off for non testing #define PRECOLLIDE_MODS(rho,ux,uy,uz) -#else // LBM_INCLUDE_TESTSOLVERS!=1 - -#define PRECOLLIDE_GETPOS \ - LbmVec( \ - ((this->mvGeoEnd[0]-this->mvGeoStart[0])/(LbmFloat)this->mSizex) * (LbmFloat)i + this->mvGeoStart[0], \ - ((this->mvGeoEnd[1]-this->mvGeoStart[1])/(LbmFloat)this->mSizey) * (LbmFloat)j + this->mvGeoStart[1], \ - ((this->mvGeoEnd[2]-this->mvGeoStart[2])/(LbmFloat)this->mSizez) * (LbmFloat)k + this->mvGeoStart[2] \ - ) - -// debug modifications of collide vars (testing) -#define PRECOLLIDE_MODS(rho,ux,uy,uz) \ - if( (this->mTForceStrength>0.) && (RFLAG(0,i,j,k, mLevel[0].setCurr)&CFNoBndFluid) ) { \ - LbmVec pos = PRECOLLIDE_GETPOS; \ - LbmVec vel(ux,uy,uz); \ - mpTest->mControlParts.modifyVelocity(pos,vel); /* */\ - if((i==16)&&(j==10)){ /*debugMarkCell(0,16,10,0);*/ errMsg("FTDEB"," at "<<PRINT_IJK<<" targ:"<<vel<<",len:"<<norm(vel)<<", org:"<<ntlVec3Gfx(ux,uy,uz) ); }\ - ux = vel[0]; uy=vel[1]; uz=vel[2]; \ - /* test acutal values...? */ \ - } - /* end PRECOLLIDE_MODS */ - -// debug modifications of collide vars (testing) -#define _PRECOLLIDE_MODS(rho,ux,uy,uz) \ - if(this->mTForceStrength>0.) { \ - LbmVec u(ux,uy,uz); \ - LbmVec pos = PRECOLLIDE_GETPOS; \ - LbmVec targv = u; const int lev=0; \ - mpTest->mControlParts.modifyVelocity(pos,targv); /* */\ - LbmVec devia = (targv-u); \ - LbmVec set; \ - set = u + (devia/this->mOmega) * this->mTForceStrength; /* */\ - set = u + targv*this->mTForceStrength; /* */\ - ux = set[0]; uy=set[1]; uz=set[2]; \ - if(j==10){ errMsg("FTDEB"," at "<<PRINT_IJK<<" targ"<<targv<<","<<norm(targv)<<" set:"<<set<<","<<norm(set) ); }\ - } - /* end PRECOLLIDE_MODS */ - -#endif // LBM_INCLUDE_TESTSOLVERS!=1 - -/* - - \ - if((0) && (j>=0.01*this->mSizey) && (j<0.9*this->mSizey)) { \ - if((1) && (i>=0.5*this->mSizex) && (i<0.6*this->mSizex)) { \ - LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \ - LbmVec u(ux,uy,uz); LbmVec t(1., 0., 0.); \ - ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; \ - } \ - if((1) && (i>=0.65*this->mSizex) && (i<0.75*this->mSizex)) { \ - LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \ - LbmVec u(ux,uy,uz); LbmVec t(0., 1., 0.); \ - ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; \ - } \ - } \ - \ - if((0) && (j>=0.0*this->mSizey) && (j<1.0 *this->mSizey)) { \ - if((1) && (i>=0.6 *this->mSizex) && (i<0.7 *this->mSizex)) { \ - LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \ - LbmVec u(ux,uy,uz); \ - LbmVec t(0., 1., 0.); \ - LbmFloat devia = norm(u-t); \ - // 0.0001-strong, 0.000001-small, - t = u - (devia/this->mOmega) * this->mTForceStrength; \ - // ux = len*t[0]; uy=len*t[1]; uz=len*t[2]; - ux = t[0]; uy=t[1]; uz=t[2]; \ - //ux= uy= uz= 0.0; - //errMsg("DDDD"," at "<<PRINT_IJK); - } \ - } \ - if(0) { \ - LbmFloat len = sqrtf(ux*ux + uy*uy + uz*uz); \ - LbmVec u(ux,uy,uz); \ - LbmVec p(i/(LbmFloat)this->mSizex, j/(LbmFloat)this->mSizey, k/(LbmFloat)this->mSizez); \ - LbmVec cp(0.5, 0.5, 0.5); \ - LbmVec delt = cp - p; \ - LbmFloat dist = norm(delt); \ - normalize(delt); \ - LbmVec tang = cross(delt, LbmVec(0.,0.,1.)); \ - normalize(tang); \ - const LbmFloat falloff = 5.0; \ - LbmVec targv = tang*1.0 * 1.0/(1.0+falloff*dist); \ - LbmVec devia = (targv-u); \ - LbmFloat devial = norm(devia); \ - LbmVec set; \ - set = u +targv * this->mTForceStrength; \ - set = u + (devia/this->mOmega) * this->mTForceStrength; \ - ux = set[0]; uy=set[1]; uz=set[2]; \ - if(j==10){ errMsg("FTDEB"," at "<<PRINT_IJK<<" dist:"<<delt<<","<<dist<<" targ"<<targv<<","<<norm(targv)<<" set:"<<set<<","<<norm(set) ); }\ - } \ - -*/ /****************************************************************************** * normal relaxation @@ -1152,14 +1059,14 @@ inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmF } #define DEBUG_CALCPRINTCELL(str,df) {\ - LbmFloat rho=df[0], ux=0., uy=0., uz=0.; \ + LbmFloat prho=df[0], pux=0., puy=0., puz=0.; \ for(int l=1; l<this->cDfNum; l++) { \ - rho += df[l]; \ - ux += (this->dfDvecX[l]*df[l]); \ - uy += (this->dfDvecY[l]*df[l]); \ - uz += (this->dfDvecZ[l]*df[l]); \ + prho += df[l]; \ + pux += (this->dfDvecX[l]*df[l]); \ + puy += (this->dfDvecY[l]*df[l]); \ + puz += (this->dfDvecZ[l]*df[l]); \ } \ - errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<rho<<" vel="<<ntlVec3Gfx(ux,uy,uz) ); \ + errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<prho<<" vel="<<ntlVec3Gfx(pux,puy,puz) ); \ } /* END DEBUG_CALCPRINTCELL */ // "normal" collision diff --git a/intern/elbeem/intern/solver_util.cpp b/intern/elbeem/intern/solver_util.cpp index 7406c55acbc..8cd698eb65f 100644 --- a/intern/elbeem/intern/solver_util.cpp +++ b/intern/elbeem/intern/solver_util.cpp @@ -110,39 +110,6 @@ void LbmFsgrSolver::prepareVisualization( void ) { *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); } - -#if LBM_INCLUDE_TESTSOLVERS==1 - /*if(mUseTestdata) { - int border = 1; - for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++) - for(int j=0;j<mLevel[mMaxRefine].lSizey-1;j++) { - for(int l=0; l<=border; l++) { - *this->mpIso->lbmGetData( l-1, j,ZKOFF) = *this->mpIso->lbmGetData( border+1, j,ZKOFF); - *this->mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-l, j,ZKOFF) = *this->mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-border-1, j,ZKOFF); - } - } - - for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++) - for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) { - for(int l=0; l<=border; l++) { - *this->mpIso->lbmGetData( i, l-1, ZKOFF) = *this->mpIso->lbmGetData( i, border+1, ZKOFF); - *this->mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-l, ZKOFF) = *this->mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-border-1, ZKOFF); - } - } - - if(LBMDIM == 3) { - // only for 3D - for(int j=-1;j<mLevel[mMaxRefine].lSizey+1;j++) - for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) { - for(int l=0; l<=border; l++) { - *this->mpIso->lbmGetData( i,j,l-1 ) = *this->mpIso->lbmGetData( i,j, border+1 ); - *this->mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-l) = *this->mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-1-border); - } - } - } - } // testdata */ -#endif // LBM_INCLUDE_TESTSOLVERS==1 - // */ // update preview, remove 2d? if((this->mOutputSurfacePreview)&&(LBMDIM==3)) { @@ -175,7 +142,7 @@ void LbmFsgrSolver::prepareVisualization( void ) { *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , this->mSizez-1); } - if(mUseTestdata) { + if(mFarFieldSize>=2.) { // also remove preview border for(int k= 0; k< (pvsz-1); ++k) { for(int j=0;j< pvsy;j++) { @@ -267,11 +234,11 @@ vector<ntlGeometryObject*> LbmFsgrSolver::getDebugObjects() { *****************************************************************************/ /*! init particle positions */ -int LbmFsgrSolver::initParticles(ParticleTracer *partt) { +int LbmFsgrSolver::initParticles() { int workSet = mLevel[mMaxRefine].setCurr; int tries = 0; int num = 0; - mpParticles=partt; + ParticleTracer *partt = mpParticles; partt->setStart( this->mvGeoStart + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) ); partt->setEnd ( this->mvGeoEnd + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) ); @@ -279,11 +246,14 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) { partt->setSimStart( ntlVec3Gfx(0.0) ); partt->setSimEnd ( ntlVec3Gfx(this->mSizex, this->mSizey, getForZMaxBnd(mMaxRefine)) ); - while( (num<partt->getNumParticles()) && (tries<100*partt->getNumParticles()) ) { + while( (num<partt->getNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) { LbmFloat x,y,z; - x = 0.0+(( (LbmFloat)(this->mSizex-1) ) * (rand()/(RAND_MAX+1.0)) ); - y = 0.0+(( (LbmFloat)(this->mSizey-1) ) * (rand()/(RAND_MAX+1.0)) ); - z = 0.0+(( (LbmFloat) getForZMax1(mMaxRefine) )* (rand()/(RAND_MAX+1.0)) ); + //x = 0.0+(( (LbmFloat)(this->mSizex-1) ) * (rand()/(RAND_MAX+1.0)) ); + //y = 0.0+(( (LbmFloat)(this->mSizey-1) ) * (rand()/(RAND_MAX+1.0)) ); + //z = 0.0+(( (LbmFloat) getForZMax1(mMaxRefine) )* (rand()/(RAND_MAX+1.0)) ); + x = 1.0+(( (LbmFloat)(this->mSizex-3.) ) * (rand()/(RAND_MAX+1.0)) ); + y = 1.0+(( (LbmFloat)(this->mSizey-3.) ) * (rand()/(RAND_MAX+1.0)) ); + z = 1.0+(( (LbmFloat) getForZMax1(mMaxRefine)-2. )* (rand()/(RAND_MAX+1.0)) ); int i = (int)(x+0.5); int j = (int)(y+0.5); int k = (int)(z+0.5); @@ -292,12 +262,19 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) { z = 0.5; // place in the middle of domain } - if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) || - TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { // only fluid cells? + //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) || + //TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { // only fluid cells? + if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) + //&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid ) + ) { // inner fluid only + bool cellOk = true; + FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; } + if(!cellOk) continue; // in fluid... partt->addParticle(x,y,z); partt->getLast()->setStatus(PART_IN); - partt->getLast()->setType(PART_BUBBLE); + partt->getLast()->setType(PART_TRACER); + partt->getLast()->setSize(1.0); num++; } tries++; @@ -337,6 +314,7 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) { if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() ); } } + // place floats on rectangular region FLOAT_JITTER_BND if(mpTest->mDebugvalue2==-10.0){ const int lev = mMaxRefine; const int sx = mLevel[lev].lSizex; @@ -348,8 +326,9 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) { LbmFloat x,y,z; x = 0.0+(LbmFloat)(i); y = 0.0+(LbmFloat)(j); - //z = 0.5+(LbmFloat)(k); - z = mLevel[lev].lSizez/20.0*8.0 - 1.0; + //z = mLevel[lev].lSizez/10.0*2.5 - 1.0; + z = mLevel[lev].lSizez/20.0*7.5 - 1.0; + //z = mLevel[lev].lSizez/20.0*4.5 - 1.0; partt->addParticle(x,y,z); //if( (i>0)&&(i<sx) && (j>0)&&(j<sy) ) { partt->getLast()->setStatus(PART_IN); } else { partt->getLast()->setStatus(PART_OUT); } partt->getLast()->setStatus(PART_IN); @@ -363,17 +342,21 @@ int LbmFsgrSolver::initParticles(ParticleTracer *partt) { #endif // LBM_INCLUDE_TESTSOLVERS - debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb, 10); + debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb<<", tries:"<<tries, 10); if(num != partt->getNumParticles()) return 1; return 0; } -void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { +#define P_CHANGETYPE(p, newtype) \ + p->setLifeTime(0); \ + /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \ + p->setType(newtype); + +void LbmFsgrSolver::advanceParticles() { int workSet = mLevel[mMaxRefine].setCurr; LbmFloat vx=0.0,vy=0.0,vz=0.0; LbmFloat rho, df[27]; //feq[27]; - if(mpParticles!=partt) { errMsg("LbmFsgrSolver::advanceParticles","Invalid ParticleTracer..."); } #define DEL_PART { \ /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<" "); */ \ p->setActive( false ); \ @@ -394,6 +377,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { const LbmFloat v1 = 9.0; // v max const LbmFloat v2 = 2.0; // v min const LbmVec rwgrav = vec2L( this->mpParam->getGravity(mSimulationTime) ); + const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){ // TODO use timestep size //bool isIn,isOut,isInZ; @@ -402,11 +386,11 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { //const int cutval = 0; // TODO FIXME add half border! const int cutval = mCutoff; // use full border!? int actCnt=0; - if(this->mStepCnt%50==49) { partt->cleanup(); } - for(vector<ParticleObject>::iterator pit= partt->getParticlesBegin(); - pit!= partt->getParticlesEnd(); pit++) { - //errMsg("PIT"," pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel() ); - //errMsg("PIT"," pit pos:"<< (*pit).getPos()<<" vel:"<< (*pit).getVel()<<" status:"<<convertFlags2String((*pit).getFlags()) <<" " <<partt->getStart()<<" "<<partt->getEnd() ); + if(this->mStepCnt%50==49) { mpParticles->cleanup(); } + for(vector<ParticleObject>::iterator pit= mpParticles->getParticlesBegin(); + pit!= mpParticles->getParticlesEnd(); pit++) { + //errMsg("PIT"," pit"<<(*pit)->getId()<<" pos:"<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel() ); + //errMsg("PIT"," pit pos:"<< (*pit).getPos()<<" vel:"<< (*pit).getVel()<<" status:"<<convertFlags2String((*pit).getFlags()) <<" " <<mpParticles->getStart()<<" "<<mpParticles->getEnd() ); //int flag = (*pit).getFlags(); if( (*pit).getActive()==false ) continue; int i,j,k; @@ -418,25 +402,23 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { i= (int)(pos[0]+0.5); j= (int)(pos[1]+0.5); k= (int)(pos[2]+0.5); - if(LBMDIM==2) { - k = 0; - } + if(LBMDIM==2) { k = 0; } // only testdata handling, all for sws #if LBM_INCLUDE_TESTSOLVERS==1 - if(mUseTestdata) { - if(mpTest->mDebugvalue1>0.0){ - p->setStatus(PART_OUT); - mpTest->handleParticle(p, i,j,k); continue; - } } + if(useff) { + p->setStatus(PART_OUT); + mpTest->handleParticle(p, i,j,k); continue; + } #endif // LBM_INCLUDE_TESTSOLVERS==1 + // FIXME , add k tests again, remove per type ones... if(p->getStatus()&PART_IN) { // IN if( (i<cutval)||(i>this->mSizex-1-cutval)|| (j<cutval)||(j>this->mSizey-1-cutval) //||(k<cutval)||(k>this->mSizez-1-cutval) ) { - if(!mUseTestdata) { DEL_PART; + if(!useff) { DEL_PART; } else { p->setStatus(PART_OUT); /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART; @@ -454,7 +436,8 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } //p->setStatus(PART_OUT);// DEBUG always out! - if(p->getType()==PART_BUBBLE) { + if( (p->getType()==PART_BUBBLE) || + (p->getType()==PART_TRACER) ) { // no interpol rho = vx = vy = vz = 0.0; @@ -481,14 +464,14 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } else { // OUT // out of bounds, deactivate... // FIXME make fsgr treatment - p->setType( PART_FLOAT ); continue; + if(p->getType()==PART_BUBBLE) { P_CHANGETYPE(p, PART_FLOAT ); continue; } } } else { // below 3d region, just rise } } else { // OUT #if LBM_INCLUDE_TESTSOLVERS==1 - if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); } + if(useff) { mpTest->handleParticle(p, i,j,k); } else DEL_PART; #else // LBM_INCLUDE_TESTSOLVERS==1 DEL_PART; @@ -497,7 +480,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } ntlVec3Gfx v = p->getVel(); // dampen... - if(mUseTestdata) { + if( (useff)&& (p->getType()==PART_BUBBLE) ) { // test rise //O vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2]; @@ -533,11 +516,35 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { LbmFloat w = 0.99; vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[mMaxRefine].gravity[2]); v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2); + } else if(p->getType()==PART_TRACER) { + v = ntlVec3Gfx(vx,vy,vz); + if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid) ) { + // ok + } else { + const int lev = mMaxRefine; + LbmFloat nx,ny,nz, nv1,nv2; + //mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux); + if(RFLAG_NB(lev,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dE,dFfrac); } else nv1 = 0.0; + if(RFLAG_NB(lev,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dW,dFfrac); } else nv2 = 0.0; + nx = 0.5* (nv2-nv1); + if(RFLAG_NB(lev,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dN,dFfrac); } else nv1 = 0.0; + if(RFLAG_NB(lev,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dS,dFfrac); } else nv2 = 0.0; + ny = 0.5* (nv2-nv1); +#if LBMDIM==3 + if(RFLAG_NB(lev,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dT,dFfrac); } else nv1 = 0.0; + if(RFLAG_NB(lev,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dB,dFfrac); } else nv2 = 0.0; + nz = 0.5* (nv2-nv1); +#else //LBMDIM==3 + nz = 0.0; +#endif //LBMDIM==3 + + v = (ntlVec3Gfx(nx,ny,nz)) * -0.025; + } } + p->setVel( v ); - //p->setVel( ntlVec3Gfx(vx,vy,vz) ); p->advanceVel(); - // fluid particle + //errMsg("PPPP"," pos"<<p->getPos()<<" "<<p->getVel() ); } // drop handling @@ -583,13 +590,17 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { if(k<cutval) { DEL_PART; continue; } if(k<=this->mSizez-1-cutval){ //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter)) { - if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty)) { + if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) { // still ok - } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){ + // shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){ + } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid) ){ // FIXME make fsgr treatment - if(p->getLifeTime()>50) { - p->setType( PART_FLOAT ); continue; - } else DEL_PART; + //if(p->getLifeTime()>50) { + P_CHANGETYPE(p, PART_FLOAT ); continue; + // jitter in cell to prevent stacking when hitting a steep surface + LbmVec pos = p->getPos(); pos[0] += (rand()/(RAND_MAX+1.0))-0.5; + pos[1] += (rand()/(RAND_MAX+1.0))-0.5; p->setPos(pos); + //} else DEL_PART; } else { DEL_PART; this->mNumParticlesLost++; @@ -597,7 +608,7 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } } else { // OUT #if LBM_INCLUDE_TESTSOLVERS==1 - if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); } + if(useff) { mpTest->handleParticle(p, i,j,k); } else{ DEL_PART; } #else // LBM_INCLUDE_TESTSOLVERS==1 { DEL_PART; } @@ -619,9 +630,10 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) { //errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags()) ); - //p->setType( PART_BUBBLE ); continue; + //P_CHANGETYPE(p, PART_BUBBLE ); continue; // currently bubbles off! DEBUG! - DEL_PART; // DEBUG bubbles off for siggraph + //DEL_PART; // DEBUG bubbles off for siggraph + P_CHANGETYPE(p, PART_FLOAT ); continue; } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { //errMsg("PIT","NEWDROP pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags()) ); //? if(p->getLifeTime()>50) { @@ -631,13 +643,14 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { //(j<=cutval)||(j>=this->mSizey-1-cutval)) { //} else //if(p->getLifeTime()>10) { - p->setType( PART_DROP ); continue; + P_CHANGETYPE(p, PART_DROP ); continue; //} else DEL_PART; } } else { // OUT // undecided particle outside... remove? DEL_PART; + //P_CHANGETYPE(p, PART_FLOAT ); continue; } } @@ -650,11 +663,17 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { #if LBM_INCLUDE_TESTSOLVERS==1 // vanishing prob = (rand()/(RAND_MAX+1.0)); - if((mUseTestdata)&&(k>mpTest->mFluidHeight)) { - LbmFloat fac = (LbmFloat)(k-mpTest->mFluidHeight)/(LbmFloat)(10*(mLevel[mMaxRefine].lSizez-mpTest->mFluidHeight)); - prob /= fac; // TODO test? errMsg("T","T "<<prob<<" "<<fac); + // increse del prob. up to max height by given factor + const int fhCheckh = (int)(mpTest->mFluidHeight*1.25); + const LbmFloat fhProbh = 25.; + if((useff)&&(k>fhCheckh)) { + LbmFloat fac = (LbmFloat)(k-fhCheckh)/(LbmFloat)(fhProbh*(mLevel[mMaxRefine].lSizez-fhCheckh)); + prob /= fac; } if(prob<mLevel[mMaxRefine].timestep*0.1) DEL_PART; + // forced vanishing + //? if(k>this->mSizez*3/4) { if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;} + #else // LBM_INCLUDE_TESTSOLVERS==1 #endif // LBM_INCLUDE_TESTSOLVERS==1 @@ -663,13 +682,15 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { //ntlVec3Gfx v = getVelocityAt(i,j,k); rho = vx = vy = vz = 0.0; - //const int DEPTH_AVG=11; // TODO how much!? - const int DEPTH_AVG=7; // TODO how much!? + // define from particletracer.h +#if MOVE_FLOATS==1 + const int DEPTH_AVG=3; // only average interface vels int ccnt=0; - for(int kk=1;kk<DEPTH_AVG;kk+=2) { + for(int kk=0;kk<DEPTH_AVG;kk+=1) { //for(int kk=1;kk<DEPTH_AVG;kk+=1) { if((k-kk)<1) continue; - if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue; + //if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue; + if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue; ccnt++; FORDF0{ LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l); @@ -681,32 +702,44 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } } if(ccnt) { - vx /=(LbmFloat)(ccnt * 1.0); // half xy speed! value2 - vy /=(LbmFloat)(ccnt * 1.0); + // use halved surface velocity (todo, use omega instead) + vx /=(LbmFloat)(ccnt * 2.0); // half xy speed! value2 + vy /=(LbmFloat)(ccnt * 2.0); vz /=(LbmFloat)(ccnt); } - // forced vanishing - if(k>this->mSizez*3/4) { if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;} +#else // MOVE_FLOATS==1 + vx=vy=0.; //p->setVel(ntlVec3Gfx(0.) ); // static_float +#endif // MOVE_FLOATS==1 + vx += (rand()/(RAND_MAX+1.0))*FLOAT_JITTER-(FLOAT_JITTER*0.5); + vy += (rand()/(RAND_MAX+1.0))*FLOAT_JITTER-(FLOAT_JITTER*0.5); + bool delfloat = false; if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) { + // in fluid cell if((1) && (k<this->mSizez-3) && ( TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFInter ) || - TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) ) - ) { + TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) + ) ) { vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2]; if(vz<0.0) vz=0.0; - } else DEL_PART; + } else delfloat=true; + // keep below obstacles + if((delfloat) && (TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFBnd )) ) { + delfloat=false; vz=0.0; + } } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { // keep in interface , one grid cell offset is added in part. gen - } - // check if above inter, remove otherwise - else if((1) && (k>2) && ( - TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) || - TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) ) - ) { - vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2]; - if(vz>0.0) vz=0.0; - } else DEL_PART; // */ + } else { //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { // shipt?, was all else before + // check if above inter, remove otherwise + if((1) && (k>2) && ( + TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) || + TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) + ) ) { + vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2]; + if(vz>0.0) vz=0.0; + } else delfloat=true; // */ + } + if(delfloat) DEL_PART; /* // move down from empty else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { @@ -715,20 +748,38 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { //DEL_PART; // ???? } else { DEL_PART; } // */ //vz = 0.0; // DEBUG - ntlVec3Gfx v(vx,vy,vz); - p->setVel( vec2G(v) ); //? + + p->setVel( vec2G( ntlVec3Gfx(vx,vy,vz) ) ); //? //p->setVel( vec2G(v)*0.75 + p->getVel()*0.25 ); //? p->advanceVel(); //errMsg("PIT","IN pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags()) ); } else { #if LBM_INCLUDE_TESTSOLVERS==1 - if(mUseTestdata) { mpTest->handleParticle(p, i,j,k); } + if(useff) { mpTest->handleParticle(p, i,j,k); } else DEL_PART; #else // LBM_INCLUDE_TESTSOLVERS==1 DEL_PART; #endif // LBM_INCLUDE_TESTSOLVERS==1 //errMsg("PIT","OUT pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags()) ); } + + // additional bnd jitter + if((1) && (useff) && (p->getLifeTime()<3)) { + // use half butoff border 1/8 + int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5); + if(maxdw<3) maxdw=3; +#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0) +#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) + //if(ABS(i-( cutval))<maxdw) { p->advance( (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5) ,0.,0.); } + //if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance( (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5) ,0.,0.); } + if((j>=0)&&(j<=this->mSizey-1)) { + if(ABS(i-( cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-( cutval))), 0.,0.); } + if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); } + } +//#undef FLOAT_JITTER_BND +#undef FLOAT_JITTBNDRAND + //if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval) + } } // unknown particle type @@ -737,12 +788,13 @@ void LbmFsgrSolver::advanceParticles(ParticleTracer *partt) { } } myTime_t parttend = getTime(); - debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<" "<<partt->getNumParticles() , 10 ); + debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<" "<<mpParticles->getNumParticles() , 10 ); } -void LbmFsgrSolver::notifySolverOfDump(int frameNr,char *frameNrStr,string outfilename) { +void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) { // 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"; @@ -752,7 +804,9 @@ void LbmFsgrSolver::notifySolverOfDump(int frameNr,char *frameNrStr,string outfi 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 = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); + 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) & CFFluid) val = 1.; //fwrite( &val, sizeof(val), 1, file); // binary fprintf(file, "%f ",val); // text //errMsg("W", PRINT_IJK<<" val:"<<val); |