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
This commit is contained in:
2006-03-29 07:35:54 +00:00
parent 0d2902b1fe
commit 0a63b3c0ca
30 changed files with 622 additions and 434 deletions

View File

@@ -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

View File

@@ -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

View File

@@ -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) {

View File

@@ -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

View File

@@ -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) {

View File

@@ -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() {

View File

@@ -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

View File

@@ -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
}
/*****************************************************************************/

View File

@@ -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);

View File

@@ -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:

View File

@@ -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;
};

View File

@@ -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,14 +186,16 @@ 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();
@@ -225,6 +203,9 @@ void ParticleTracer::notifyOfDump(int frameNr,char *frameNrStr,string outfilenam
// 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

View File

@@ -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;

View File

@@ -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);
}

View File

@@ -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();

View File

@@ -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)])

View File

@@ -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;
}

View File

@@ -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);
}

View File

@@ -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;
};

View File

@@ -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--;
}
}

View File

@@ -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

View File

@@ -111,39 +111,6 @@ void LbmFsgrSolver::prepareVisualization( void ) {
}
#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)) {
int pvsx = (int)(this->mPreviewFactor*this->mSizex);
@@ -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
} 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; // */
}
// 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; // */
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);

View File

@@ -2392,7 +2392,7 @@ void readVelgz(char *filename, Object *srcob)
vverts[i].co[j] = 0.;
}
}
if(srcob->fluidsimSettings->typeFlags&OB_FSDOMAIN_NOVECGEN) return;
if(srcob->fluidsimSettings->domainNovecgen>0) return;
if(len<7) {
//printf("readVelgz Eror: invalid filename '%s'\n",filename); // DEBUG
@@ -2516,6 +2516,18 @@ void loadFluidsimMesh(Object *srcob, int useRenderParams)
mesh = readBobjgz(targetFile, srcob->fluidsimSettings->orgMesh, bbSize,bbSize );
}
if(!mesh) {
// switch, abort background rendering when fluidsim mesh is missing
const char *strEnvName2 = "BLENDER_ELBEEMBOBJABORT"; // from blendercall.cpp
if(G.background==1) {
if(getenv(strEnvName2)) {
int elevel = atoi(getenv(strEnvName2));
if(elevel>0) {
printf("Env. var %s set, fluid sim mesh not found, aborting render...\n",strEnvName2);
exit(1);
}
}
}
// display org. object upon failure
srcob->data = srcob->fluidsimSettings->orgMesh;
return;

View File

@@ -1578,6 +1578,11 @@ static pMatrixCache *cache_object_matrices(Object *ob, int start, int end)
return mcache;
}
/* for fluidsim win32 debug messages */
#if defined(WIN32) && (!(defined snprintf))
#define snprintf _snprintf
#endif
/* main particle building function
one day particles should become dynamic (realtime) with the current method as a 'bake' (ton) */
void build_particle_system(Object *ob)
@@ -1599,7 +1604,7 @@ void build_particle_system(Object *ob)
float *volengths= NULL, *folengths= NULL;
int deform=0, a, totpart, paf_sta, paf_end;
int waitcursor_set= 0, totvert, totface, curface, curvert;
int readMask =0, activeParts;
int readMask, activeParts, fileParts;
/* return conditions */
if(ob->type!=OB_MESH) return;
@@ -1617,6 +1622,7 @@ void build_particle_system(Object *ob)
char *suffix = "fluidsurface_particles_#";
char *suffix2 = ".gz";
char filename[256];
char debugStrBuffer[256];
int curFrame = G.scene->r.cfra -1; // warning - sync with derived mesh fsmesh loading
int j, numFileParts;
gzFile gzf;
@@ -1635,8 +1641,7 @@ void build_particle_system(Object *ob)
gzf = gzopen(filename, "rb");
if (!gzf) {
//char debugStrBuffer[256];
//define win32... snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s'\n", filename);
snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename);
//elbeemDebugOut(debugStrBuffer);
paf->totpart = 1;
return;
@@ -1648,28 +1653,24 @@ void build_particle_system(Object *ob)
paf->totpart= totpart;
paf->totkey= 1;
/* initialize particles */
new_particle(paf);// ?
ftime = 0.0;
dtime= 0.0f;
new_particle(paf);
ftime = 0.0; // unused...
// set up reading mask
//for(j=1; j<=4; j++ ){ if(ob->fluidsimSettings->guiDisplayMode&j) readMask |= (1<<j); }
readMask = ob->fluidsimSettings->guiDisplayMode;
readMask = ob->fluidsimSettings->typeFlags;
activeParts=0;
// FIXME only allocate needed ones?
fileParts=0;
//fprintf(stderr,"FSPARTICLE debug set %s , tot%d mask=%d \n", filename, totpart, readMask );
for(a=0; a<totpart; a++, ftime+=dtime) {
for(a=0; a<totpart; a++) {
int ptype=0;
short shsize=0;
float convertSize=0.0;
gzread(gzf, &ptype, sizeof( ptype ));
//if(a<25) fprintf(stderr,"FSPARTICLE debug set %s , a%d t=%d , mask=%d , active%d\n", filename, a, ptype, readMask, activeParts );
if(ptype&readMask) {
activeParts++;
pa= new_particle(paf);
pa->time= ftime;
pa->lifetime= ftime + G.scene->r.efra +1.0;
pa->lifetime= ftime + 10000.; // add large number to make sure they are displayed, G.scene->r.efra +1.0;
pa->co[0] = 0.0;
pa->co[1] =
pa->co[2] = 1.0*(float)a / (float)totpart;
@@ -1692,18 +1693,20 @@ void build_particle_system(Object *ob)
gzread(gzf, &wrf, sizeof( wrf ));
vel[j] = wrf;
}
//if(a<25) fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f \n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
} else {
// skip...
for(j=0; j<2*3+1; j++) {
float wrf; gzread(gzf, &wrf, sizeof( wrf ));
}
}
//fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f \n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
fileParts++;
}
gzclose( gzf );
totpart = paf->totpart = activeParts;
//fprintf(stderr,"PARTOBH debug %s %d \n", ob->id.name, totpart); // DEBUG
snprintf(debugStrBuffer,256,"readFsPartData::done - particles:%d, active:%d, file:%d, mask:%d \n", paf->totpart,activeParts,fileParts,readMask);
elbeemDebugOut(debugStrBuffer);
return;
}

View File

@@ -2531,6 +2531,7 @@ static void direct_link_object(FileData *fd, Object *ob)
ob->fluidsimSettings->orgMesh = NULL; //ob->data;
ob->fluidsimSettings->meshSurface = NULL;
ob->fluidsimSettings->meshBB = NULL;
ob->fluidsimSettings->meshSurfNormals = NULL;
}
link_list(fd, &ob->prop);

View File

@@ -96,14 +96,19 @@ typedef struct FluidsimSettings {
/* additional flags depending on the type, lower short contains flags
* to check validity, higher short additional flags */
int typeFlags;
short typeFlags;
char domainNovecgen,dummyc;
/* boundary "stickiness" for part slip values */
float partSlipValue;
/* particle generation - on if >0, then determines amount */
float generateParticles, dummy;
float generateParticles;
/* smooth fluid surface? */
float surfaceSmoothing;
/* particle display - size scaling, and alpha influence */
float particleInfSize, particleInfAlpha;
/* testing vars */
float farFieldSize,dummyf;
/* save fluidsurface normals in mvert.no, and surface vertex velocities (if available) in mvert.co */
struct MVert *meshSurfNormals;
@@ -119,13 +124,12 @@ typedef struct FluidsimSettings {
#define OB_FLUIDSIM_OUTFLOW 32
#define OB_FLUIDSIM_PARTICLE 64
#define OB_TYPEFLAG_START 16
#define OB_TYPEFLAG_START 0
#define OB_FSGEO_THIN (1<<(OB_TYPEFLAG_START+1))
#define OB_FSBND_NOSLIP (1<<(OB_TYPEFLAG_START+2))
#define OB_FSBND_PARTSLIP (1<<(OB_TYPEFLAG_START+3))
#define OB_FSBND_FREESLIP (1<<(OB_TYPEFLAG_START+4))
#define OB_FSINFLOW_LOCALCOORD (1<<(OB_TYPEFLAG_START+5))
#define OB_FSDOMAIN_NOVECGEN (1<<(OB_TYPEFLAG_START+6))
// guiDisplayMode particle flags
#define OB_FSDOM_GEOM 1

View File

@@ -927,7 +927,7 @@ static void render_particle_system(Render *re, Object *ob, PartEff *paf)
if(useFluidsimParticles) {
// rescale to 1.0-10.0, then div by 5 afterwards, gives values in range 0.2-2.0
double fspsize = ((double)pa->rt / 1000.0f) / 5.0 ;
haloScale = (float)pow(fspsize, (double)ob->fluidsimSettings->particleInfSize);
haloScale = 1.0/(float)pow(fspsize, (double)ob->fluidsimSettings->particleInfSize);
ma->alpha = iniAlpha / (float)pow( fspsize, (double)ob->fluidsimSettings->particleInfAlpha);
if(ma->alpha>1.) ma->alpha = 1.;
}

View File

@@ -2454,25 +2454,31 @@ static void object_panel_fluidsim(Object *ob)
uiDefBut(block, LABEL, 0, "Domain boundary type settings:", 0,yline,300,objHeight, NULL, 0.0, 0, 0, 0, "");
yline -= lineHeight;
uiBlockBeginAlign(block);
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Noslip", 00, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_NOSLIP, 20.0, 1.0, "Obstacle causes zero normal and tangential velocity (=sticky). Default for all. Only option for moving objects.");
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Part", 100, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_PARTSLIP, 20.0, 2.0, "Mix between no-slip and free-slip. Non moving objects only!");
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Free", 200, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_FREESLIP, 20.0, 3.0, "Obstacle only causes zero normal velocity (=not sticky). Non moving objects only!");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Noslip", 0, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_NOSLIP, 20.0, 1.0, "Obstacle causes zero normal and tangential velocity (=sticky). Default for all. Only option for moving objects.");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Part", 100, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_PARTSLIP, 20.0, 2.0, "Mix between no-slip and free-slip. Non moving objects only!");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Free", 200, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_FREESLIP, 20.0, 3.0, "Obstacle only causes zero normal velocity (=not sticky). Non moving objects only!");
uiBlockEndAlign(block);
yline -= lineHeight;
uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, "");
if(fss->typeFlags&OB_FSBND_PARTSLIP) {
uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButF(block, NUM, B_DIFF, "", 150, yline,150,objHeight, &fss->partSlipValue, 0.0, 0.1, 10,0, ".");
yline -= lineHeight;
}
} else { uiDefBut(block, LABEL, 0, "-", 150,yline,150,objHeight, NULL, 0.0, 0, 0, 0, ""); }
yline -= lineHeight;
// copied from obstacle...
uiDefBut(block, LABEL, 0, "Generate Particles:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButF(block, NUM, B_DIFF, "", 200, yline,100,objHeight, &fss->generateParticles, 0.0, 10.0, 10,0, "Amount of particles to generate (0=off, 1=normal, >1=more).");
yline -= lineHeight;
uiDefBut(block, LABEL, 0, "Generate&Use SpeedVecs:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButBitI(block, TOG, OB_FSDOMAIN_NOVECGEN, REDRAWBUTSOBJECT, "Disable", 200, yline,100,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Default is to generate and use fluidsim vertex speed vectors, this option switches calculation off during bake, and disables loading.");
uiDefBut(block, LABEL, 0, "Surface Smoothing:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButF(block, NUM, B_DIFF, "", 200, yline,100,objHeight, &fss->surfaceSmoothing, 0.0, 5.0, 10,0, "Amount of surface smoothing (0=off, 1=normal, >1=stronger smoothing).");
yline -= lineHeight;
// use new variable...
uiDefBut(block, LABEL, 0, "Generate&Use SpeedVecs:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButBitC(block, TOG, 1, REDRAWBUTSOBJECT, "Disable", 200, yline,100,objHeight, &fss->domainNovecgen, 0, 0, 0, 0, "Default is to generate and use fluidsim vertex speed vectors, this option switches calculation off during bake, and disables loading.");
yline -= lineHeight;
// copied from obstacle...
}
}
else if(
@@ -2492,7 +2498,7 @@ static void object_panel_fluidsim(Object *ob)
if(fss->type == OB_FLUIDSIM_INFLOW) {
uiDefBut(block, LABEL, 0, "Local Inflow Coords", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButBitI(block, TOG, OB_FSINFLOW_LOCALCOORD, REDRAWBUTSOBJECT, "Enable", 200, yline,100,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Use local coordinates for inflow (e.g. for rotating objects).");
uiDefButBitS(block, TOG, OB_FSINFLOW_LOCALCOORD, REDRAWBUTSOBJECT, "Enable", 200, yline,100,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Use local coordinates for inflow (e.g. for rotating objects).");
yline -= lineHeight;
}
}
@@ -2504,31 +2510,26 @@ static void object_panel_fluidsim(Object *ob)
yline -= lineHeight + 5;
uiBlockBeginAlign(block);
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Noslip", 00, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_NOSLIP, 20.0, 1.0, "Obstacle causes zero normal and tangential velocity (=sticky). Default for all. Only option for moving objects.");
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Part", 100, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_PARTSLIP, 20.0, 2.0, "Mix between no-slip and free-slip. Non moving objects only!");
uiDefButI(block, ROW, REDRAWBUTSOBJECT ,"Free", 200, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_FREESLIP, 20.0, 3.0, "Obstacle only causes zero normal velocity (=not sticky). Non moving objects only!");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Noslip", 00, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_NOSLIP, 20.0, 1.0, "Obstacle causes zero normal and tangential velocity (=sticky). Default for all. Only option for moving objects.");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Part", 100, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_PARTSLIP, 20.0, 2.0, "Mix between no-slip and free-slip. Non moving objects only!");
uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Free", 200, yline,100,objHeight, &fss->typeFlags, 15.0, OB_FSBND_FREESLIP, 20.0, 3.0, "Obstacle only causes zero normal velocity (=not sticky). Non moving objects only!");
uiBlockEndAlign(block);
yline -= lineHeight;
uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, "");
if(fss->typeFlags&OB_FSBND_PARTSLIP) {
uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefButF(block, NUM, B_DIFF, "", 150, yline,150,objHeight, &fss->partSlipValue, 0.0, 0.1, 10,0, ".");
yline -= lineHeight;
}
} else { uiDefBut(block, LABEL, 0, "-", 150,yline,150,objHeight, NULL, 0.0, 0, 0, 0, ""); }
yline -= lineHeight;
yline -= lineHeight;
}
else if(fss->type == OB_FLUIDSIM_PARTICLE) {
if(fss->guiDisplayMode==0) fss->guiDisplayMode=2; // default drops
uiDefBut(block, LABEL, 0, "Part.-Type:", 0,yline,100,objHeight, NULL, 0.0, 0, 0, 0, "");
// TODO make toggle buttons
//uiDefButS(block, MENU, B_FLUIDSIM_FORCEREDRAW, "Gui%t|Bubble %x2|Drop %x4|Newparts %x8|Float %x16",
//100,yline,200,objHeight, &fss->guiDisplayMode, 0, 0, 0, 0, "Which type of particles to display.");
//uiDefButS(block, MENU, B_DIFF, "Render%t|Geometry %x1|Preview %x2|Final %x3",
//195,yline,105,objHeight, &fss->renderDisplayMode, 0, 0, 0, 0, "How to display the fluid mesh for rendering.");
uiDefBut(block, LABEL, 0, "Drops", 100,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
fss->guiDisplayMode = 4; // fix to drops for now
if(fss->typeFlags==0) fss->typeFlags=4; // default drops
fss->typeFlags = 4; // fix to drops for now
uiDefBut(block, LABEL, 0, "Part.-Type:", 0,yline,100,objHeight, NULL, 0.0, 0, 0, 0, "");
uiDefBut(block, LABEL, 0, "Drops", 100,yline,200,objHeight, NULL, 0.0, 0, 0, 0, "");
yline -= lineHeight;
uiDefBut(block, LABEL, 0, "Size Influence:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, "");

View File

@@ -84,6 +84,7 @@
#include "BSE_headerbuttons.h"
#include "mydevice.h"
#include "blendef.h"
#include "SDL.h"
#include "SDL_thread.h"
#include "SDL_mutex.h"
@@ -109,9 +110,6 @@ void initElbeemMesh(struct Object *ob, int *numVertices, float **vertices, int *
extern int start_progress_bar(void);
extern void end_progress_bar(void);
extern int progress_bar(float done, char *busy_info);
// global solver state
extern int gElbeemState;
extern char gElbeemErrorString[];
double fluidsimViscosityPreset[6] = {
-1.0, /* unused */
@@ -197,9 +195,11 @@ FluidsimSettings *fluidsimSettingsNew(struct Object *srcob)
fluidsimGetAxisAlignedBB(srcob->data, srcob->obmat, fss->bbStart, fss->bbSize, &fss->meshBB);
fss->typeFlags = 0;
fss->domainNovecgen = 0;
fss->partSlipValue = 0.0;
fss->generateParticles = 0.0;
fss->surfaceSmoothing = 1.0;
fss->particleInfSize = 0.0;
fss->particleInfAlpha = 0.0;
@@ -445,6 +445,24 @@ void fluidsimBake(struct Object *ob)
return;
}
/* no object pointer, find in selected ones.. */
if(!ob) {
Base *base;
for(base=G.scene->base.first; base; base= base->next) {
if ( ((base)->flag & SELECT)
// ignore layer setting for now? && ((base)->lay & G.vd->lay)
) {
if((!ob)&&(base->object->fluidsimFlag & OB_FLUIDSIM_ENABLE)&&(base->object->type==OB_MESH)) {
if(base->object->fluidsimSettings->type == OB_FLUIDSIM_DOMAIN) {
ob = base->object;
}
}
}
}
// no domains found?
if(!ob) return;
}
/* check if there's another domain... */
for(obit= G.main->object.first; obit; obit= obit->id.next) {
if((obit->fluidsimFlag & OB_FLUIDSIM_ENABLE)&&(obit->type==OB_MESH)) {
@@ -764,6 +782,8 @@ void fluidsimBake(struct Object *ob)
fsset.gstar = domainSettings->gstar;
fsset.maxRefine = domainSettings->maxRefine; // check <-> gridlevels
fsset.generateParticles = domainSettings->generateParticles;
fsset.surfaceSmoothing = domainSettings->surfaceSmoothing;
fsset.farFieldSize = domainSettings->farFieldSize;
strcpy( fsset.outputPath, targetFile);
// domain channels
@@ -778,8 +798,7 @@ void fluidsimBake(struct Object *ob)
else if((domainSettings->typeFlags&OB_FSBND_PARTSLIP)) fsset.obstacleType = FLUIDSIM_OBSTACLE_PARTSLIP;
else if((domainSettings->typeFlags&OB_FSBND_FREESLIP)) fsset.obstacleType = FLUIDSIM_OBSTACLE_FREESLIP;
fsset.obstaclePartslip = domainSettings->partSlipValue;
fsset.generateVertexVectors = (int)(!(domainSettings->typeFlags&OB_FSDOMAIN_NOVECGEN));
// fprintf(stderr," VVV %d %d \n",fsset.generateVertexVectors , (domainSettings->typeFlags&OB_FSDOMAIN_NOVECGEN)); // DEBUG
fsset.generateVertexVectors = (domainSettings->domainNovecgen==0);
// init blender trafo matrix
// fprintf(stderr,"elbeemInit - mpTrafo:\n");
@@ -950,8 +969,8 @@ void fluidsimBake(struct Object *ob)
" size = " "%d" /* gridSize*/ "; \n"
" surfacepreview = " "%d" /* previewSize*/ "; \n"
" dump_velocities = " "%d" /* vector dump */ "; \n"
" smoothsurface = 1.0; \n"
" smoothnormals = 1.0; \n"
" smoothsurface = %f; \n" /* smoothing */
" smoothnormals = %f; \n"
" geoinitid = 1; \n" "\n"
" isovalue = 0.4900; \n"
" isoweightmethod = 1; \n" "\n" ;
@@ -959,7 +978,7 @@ void fluidsimBake(struct Object *ob)
fprintf(fileCfg, simString,
(double)domainSettings->realsize, (double)domainSettings->animStart, (double)domainSettings->gstar,
gridlevels, (int)domainSettings->resolutionxyz, (int)domainSettings->previewresxyz ,
(int)(!(domainSettings->typeFlags&OB_FSDOMAIN_NOVECGEN))
(int)(domainSettings->domainNovecgen==0), domainSettings->surfaceSmoothing
);
if((domainSettings->typeFlags&OB_FSBND_NOSLIP)) bi=0;
@@ -972,6 +991,7 @@ void fluidsimBake(struct Object *ob)
fluidsimPrintChannel(fileCfg, channelDomainViscosity,allchannelSize,"p_viscosity",CHANNEL_FLOAT);
fluidsimPrintChannel(fileCfg, channelDomainGravity, allchannelSize,"p_gravity",CHANNEL_VEC);
fprintf(fileCfg, " partgenprob = %f; \n", domainSettings->generateParticles); // debug test
fprintf(fileCfg, "\n} \n" );
}
@@ -1199,17 +1219,16 @@ void fluidsimBake(struct Object *ob)
if(!simAborted) {
char fsmessage[512];
char elbeemerr[256];
strcpy(fsmessage,"Fluidsim Bake Error: ");
// check if some error occurred
if(globalBakeState==-2) {
strcat(fsmessage,"Failed to initialize [Msg: ");
#ifndef WIN32
// msvc seems to have problem accessing the gElbeemErrorString var
strcat(fsmessage,"[Msg: ");
strcat(fsmessage,gElbeemErrorString);
strcat(fsmessage,"]");
#endif // WIN32
strcat(fsmessage,"|OK%x0");
elbeemGetErrorString(elbeemerr);
strcat(fsmessage,elbeemerr);
strcat(fsmessage,"] |OK%x0");
pupmenu(fsmessage);
} // init error
}

View File

@@ -1216,7 +1216,10 @@ static void winqreadview3dspace(ScrArea *sa, void *spacedata, BWinEvent *evt)
else if(G.qual==LR_CTRLKEY) {
if(okee("Bake all selected")) {
extern void softbody_bake(Object *ob);
extern void fluidsimBake(Object *ob);
softbody_bake(NULL);
// also bake first domain of selected objects...
fluidsimBake(NULL);
}
}
else if(G.qual==0)