diff --git a/intern/elbeem/COPYING b/intern/elbeem/COPYING index d876362de6c..2600c731161 100644 --- a/intern/elbeem/COPYING +++ b/intern/elbeem/COPYING @@ -1,7 +1,5 @@ - All code distributed as part of El'Bemm is covered by the following - version of the GNU General Public License, expcept for excerpts of - the trimesh2 package in src/isosurface.cpp (see COPYING_trimesh2 - for details). + All code distributed as part of El'Beem is covered by the following + version of the GNU General Public License. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License diff --git a/intern/elbeem/extern/elbeem.h b/intern/elbeem/extern/elbeem.h index ea21e807a6f..b3feda8bbe8 100644 --- a/intern/elbeem/extern/elbeem.h +++ b/intern/elbeem/extern/elbeem.h @@ -3,7 +3,7 @@ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method * All code distributed as part of El'Beem is covered by the version 2 of the * GNU General Public License. See the file COPYING for details. - * Copyright 2003-2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * API header */ @@ -72,13 +72,14 @@ typedef struct elbeemSimulationSettings { float *channelGravity; // vector /* boundary types and settings for domain walls */ - short obstacleType; - float obstaclePartslip; + short domainobsType; + float domainobsPartslip; /* generate speed vectors for vertices (e.g. for image based motion blur)*/ short generateVertexVectors; /* strength of surface smoothing */ float surfaceSmoothing; - // TODO add surf gen flags + /* no. of surface subdivisions */ + int surfaceSubdivs; /* global transformation to apply to fluidsim mesh */ float surfaceTrafo[4*4]; @@ -147,6 +148,8 @@ typedef struct elbeemMesh { /* boundary types and settings */ short obstacleType; float obstaclePartslip; + /* amount of force transfer from fluid to obj, 0=off, 1=normal */ + float obstacleImpactFactor; /* init volume, shell or both? use OB_VOLUMEINIT_xxx defines above */ short volumeInitType; @@ -230,5 +233,8 @@ double elbeemEstimateMemreq(int res, #define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6)) #define FGI_MBNDOUTFLOW (1<<(FGI_FLAGSTART+ 7)) +// all boundary types at once +#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW ) + #endif // ELBEEM_API_H diff --git a/intern/elbeem/intern/attributes.cpp b/intern/elbeem/intern/attributes.cpp index b7f194582d1..0cf5315161a 100644 --- a/intern/elbeem/intern/attributes.cpp +++ b/intern/elbeem/intern/attributes.cpp @@ -1,290 +1,44 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * - * configuration attribute storage class and attribute class + * DEPRECATED - replaced by elbeem API, only channels are still used * *****************************************************************************/ #include "attributes.h" #include "ntl_matrices.h" #include "elbeem.h" -#include -//! output attribute values? on=1/off=0 -#define DEBUG_ATTRIBUTES 0 - -//! output channel values? on=1/off=0 -#define DEBUG_CHANNELS 0 - /****************************************************************************** * attribute conversion functions *****************************************************************************/ bool Attribute::initChannel(int elemSize) { - if(!mIsChannel) return true; - if(mChannelInited==elemSize) { - // already inited... ok - return true; - } else { - // sanity check - if(mChannelInited>0) { - errMsg("Attribute::initChannel","Duplicate channel init!? ("< newvalue; - for(int i=0; i elem(elemSize); - for(int j=0; jmTimes[i+1]) { - errMsg("Attribute::initChannel","Invalid time at entry "< Attribute::getChannelFloat() { - vector timeVec; - vector valVec; - - if((!initChannel(1)) || (!mIsChannel)) { - timeVec.push_back( 0.0 ); - valVec.push_back( getAsFloat() ); - } else { - for(size_t i=0; i(valVec,timeVec); + return AnimChannel(); } - -//! get channel as integer value AnimChannel Attribute::getChannelInt() { - vector timeVec; - vector valVec; - - if((!initChannel(1)) || (!mIsChannel)) { - timeVec.push_back( 0.0 ); - valVec.push_back( getAsInt() ); - } else { - for(size_t i=0; i(valVec,timeVec); + return AnimChannel(); } - -//! get channel as integer value AnimChannel Attribute::getChannelVec3d() { - vector timeVec; - vector valVec; - - if((!initChannel(3)) || (!mIsChannel)) { - timeVec.push_back( 0.0 ); - valVec.push_back( getAsVec3d() ); - } else { - for(size_t i=0; i(valVec,timeVec); + return AnimChannel(); } - -//! get channel as float vector set AnimChannel Attribute::getChannelSetVec3f() { - vector timeVec; - ntlSetVec3f setv; - vector< ntlSetVec3f > valVec; - - if((!initChannel(3)) || (!mIsChannel)) { - timeVec.push_back( 0.0 ); - setv.mVerts.push_back( vec2F( getAsVec3d() ) ); - valVec.push_back( setv ); - } else { - for(size_t i=0; i(valVec,timeVec); + return AnimChannel(); } /****************************************************************************** * check if there were unknown params *****************************************************************************/ -bool AttributeList::checkUnusedParams() -{ - bool found = false; - for(map::iterator i=mAttrs.begin(); - i != mAttrs.end(); i++) { - if((*i).second) { - if(!(*i).second->getUsed()) { - errMsg("AttributeList::checkUnusedParams", "List "<::iterator i=mAttrs.begin(); - i != mAttrs.end(); i++) { - if((*i).second) { - (*i).second->setUsed(true); - } - } } /****************************************************************************** * Attribute list read functions *****************************************************************************/ int AttributeList::readInt(string name, int defaultValue, string source,string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readInt","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - return find(name)->getAsInt(); + return defaultValue; } bool AttributeList::readBool(string name, bool defaultValue, string source,string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readBool","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - return find(name)->getAsBool(); + return defaultValue; } double AttributeList::readFloat(string name, double defaultValue, string source,string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readFloat","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - return find(name)->getAsFloat(); + return defaultValue; } string AttributeList::readString(string name, string defaultValue, string source,string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readInt","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - return find(name)->getAsString(false); + return defaultValue; } ntlVec3d AttributeList::readVec3d(string name, ntlVec3d defaultValue, string source,string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readInt","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - return find(name)->getAsVec3d(); + return defaultValue; } void AttributeList::readMat4Gfx(string name, ntlMat4Gfx defaultValue, string source,string target, bool needed, ntlMat4Gfx *mat) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readInt","Required attribute '"<getCompleteString() <<"' as type int " , 3); } - find(name)->setUsed(true); - find(name)->getAsMat4Gfx( mat ); - return; } // set that a parameter can be given, and will be ignored... bool AttributeList::ignoreParameter(string name, string source) { - if(!exists(name)) return false; - find(name)->setUsed(true); - if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Param '"<< name <<"' set but ignored... " , 3); } - return true; + return false; } // read channels AnimChannel AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<(defaultValue); } - AnimChannel ret = find(name)->getChannelInt(); - find(name)->setUsed(true); - channelSimplifyi(ret); - return ret; + return AnimChannel(defaultValue); } AnimChannel AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<(defaultValue); } - AnimChannel ret = find(name)->getChannelFloat(); - find(name)->setUsed(true); - channelSimplifyd(ret); - return ret; + return AnimChannel(defaultValue); } AnimChannel AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<(defaultValue); } - AnimChannel ret = find(name)->getChannelVec3d(); - find(name)->setUsed(true); - channelSimplifyVd(ret); - return ret; + return AnimChannel(defaultValue); } AnimChannel AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<(defaultValue); } - AnimChannel ret = find(name)->getChannelSetVec3f(); - find(name)->setUsed(true); - //channelSimplifyVf(ret); - return ret; + return AnimChannel(defaultValue); } AnimChannel AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<(defaultValue); } - AnimChannel convert = find(name)->getChannelFloat(); - find(name)->setUsed(true); - channelSimplifyd(convert); - // convert to float vals - vector vals; - for(size_t i=0; i times = convert.accessTimes(); - AnimChannel ret(vals, times); - return ret; + return AnimChannel(defaultValue); } AnimChannel AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) { - if(!exists(name)) { - if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<(defaultValue); } - - AnimChannel convert = find(name)->getChannelVec3d(); - // convert to float - vector vals; - for(size_t i=0; i times = convert.accessTimes(); - AnimChannel ret(vals, times); - find(name)->setUsed(true); - channelSimplifyVf(ret); - return ret; + return AnimChannel(defaultValue); } /****************************************************************************** * destructor *****************************************************************************/ AttributeList::~AttributeList() { - for(map::iterator i=mAttrs.begin(); - i != mAttrs.end(); i++) { - if((*i).second) { - delete (*i).second; - (*i).second = NULL; - } - } }; @@ -560,56 +128,18 @@ AttributeList::~AttributeList() { *****************************************************************************/ //! debug function, prints value -void Attribute::print() -{ - std::ostringstream ostr; - ostr << "ATTR "<< mName <<"= "; - for(size_t i=0;i0) { - for(size_t i=0;i::iterator i=mAttrs.begin(); - i != mAttrs.end(); i++) { - if((*i).second) { - (*i).second->print(); - } - } +void AttributeList::print() { } /****************************************************************************** * import attributes from other attribute list *****************************************************************************/ -void AttributeList::import(AttributeList *oal) -{ - for(map::iterator i=oal->mAttrs.begin(); - i !=oal->mAttrs.end(); i++) { - // FIXME - check freeing of copyied attributes - if((*i).second) { - Attribute *newAttr = new Attribute( *(*i).second ); - mAttrs[ (*i).first ] = newAttr; - } - } +void AttributeList::import(AttributeList *oal) { } @@ -672,7 +202,6 @@ static bool channelSimplifyScalarT(AnimChannel &channel) { int size = channel.getSize(); if(size<=1) return false; float *nchannel = new float[2*size]; - if(DEBUG_CHANNELS) errMsg("channelSimplifyf","S" << channel.printChannel() ); // convert to array for(size_t i=0; i &channel) { times.push_back( (double)(nchannel[i*2 + 1]) ); } channel = AnimChannel(vals, times); - if(DEBUG_CHANNELS) errMsg("channelSimplifyf","C" << channel.printChannel() ); } delete [] nchannel; return ret; @@ -700,7 +228,6 @@ static bool channelSimplifyVecT(AnimChannel &channel) { int size = channel.getSize(); if(size<=1) return false; float *nchannel = new float[4*size]; - if(DEBUG_CHANNELS) errMsg("channelSimplifyf","S" << channel.printChannel() ); // convert to array for(size_t i=0; i &channel) { times.push_back( (double)(nchannel[i*4 + 3]) ); } channel = AnimChannel(vals, times); - if(DEBUG_CHANNELS) errMsg("channelSimplifyf","C" << channel.printChannel() ); } delete [] nchannel; return ret; @@ -742,13 +268,7 @@ string AnimChannel::printChannel() { return ostr.str(); } // */ -//! debug function, prints to stdout if DEBUG_CHANNELS flag is enabled, used in constructors -template -void AnimChannel::debugPrintChannel() { - if(DEBUG_CHANNELS) { errMsg("channelCons"," : " << this->printChannel() ); } -} - - +// is now in header file: debugPrintChannel() // hack to force instantiation void __forceAnimChannelInstantiation() { AnimChannel< float > tmp1; diff --git a/intern/elbeem/intern/attributes.h b/intern/elbeem/intern/attributes.h index 20e5341d5fc..19355095351 100644 --- a/intern/elbeem/intern/attributes.h +++ b/intern/elbeem/intern/attributes.h @@ -1,9 +1,9 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * - * configuration attribute storage class and attribute class + * DEPRECATED - replaced by elbeem API, only channels are still used * *****************************************************************************/ @@ -38,7 +38,7 @@ class AnimChannel ~AnimChannel() { }; // get interpolated value at time t - Scalar get(double t) { + Scalar get(double t) const { if(!mInited) { Scalar null; null=(Scalar)(0.0); return null; } if(t<=mTimes[0]) { return mValue[0]; } if(t>=mTimes[mTimes.size()-1]) { return mValue[mTimes.size()-1]; } @@ -63,7 +63,7 @@ class AnimChannel }; // get uninterpolated value at time t - Scalar getConstant(double t) { + Scalar getConstant(double t) const { //errMsg("DEBB","getc"< &accessValues() { return mValue; } //! raw access of time vector @@ -110,91 +110,6 @@ class AnimChannel }; -//! A single attribute -class Attribute -{ - public: - //! Standard constructor - Attribute(string mn, vector &value, int setline,bool channel) : - mName(mn), mValue(value), - mLine(setline), mUsed(false), mIsChannel(channel), - mChannelInited(-1) { }; - //! Copy constructor - Attribute(Attribute &a) : - mName(a.mName), mValue(a.mValue), - mLine(a.mLine), mUsed(false), mIsChannel(a.mIsChannel), - mChannelInited(a.mChannelInited), - mChannel(a.mChannel), mTimes(a.mTimes) { }; - //! Destructor - ~Attribute() { /* empty */ }; - - //! set used flag - void setUsed(bool set){ mUsed = set; } - //! get used flag - bool getUsed() { return mUsed; } - - //! set channel flag - void setIsChannel(bool set){ mIsChannel = set; } - //! get channel flag - bool getIsChannel() { return mIsChannel; } - - //! get value as string - string getAsString(bool debug=false); - //! get value as integer value - int getAsInt(); - //! get value as boolean - bool getAsBool(); - //! get value as double value - double getAsFloat(); - //! get value as 3d vector - ntlVec3d getAsVec3d(); - //! get value as 4x4 matrix - void getAsMat4Gfx(ntlMatrix4x4 *mat); - - //! get channel as integer value - AnimChannel getChannelInt(); - //! get channel as double value - AnimChannel getChannelFloat(); - //! get channel as double vector - AnimChannel getChannelVec3d(); - //! get channel as float vector set - AnimChannel getChannelSetVec3f(); - - //! get the concatenated string of all value string - string getCompleteString(); - - //! debug function, prints value - void print(); - - protected: - - /*! internal - init channel before access */ - bool initChannel(int elemSize); - - /*! the attr name */ - string mName; - - /*! the attr value */ - vector mValue; - - /*! line where the value was defined in the config file (for error messages) */ - int mLine; - - /*! was this attribute used? */ - bool mUsed; - - /*! does this attribute have a channel? */ - bool mIsChannel; - /*! does this attribute have a channel? */ - int mChannelInited; - - /*! channel attribute values (first one equals mValue) */ - vector< vector< string > > mChannel; - /*! channel attr times */ - vector< double > mTimes; -}; - - // helper class (not templated) for animated meshes class ntlSetVec3f { public: @@ -211,77 +126,68 @@ class ntlSetVec3f { vector mVerts; }; + +// warning: DEPRECATED - replaced by elbeem API +class Attribute +{ + public: + Attribute(string mn, vector &value, int setline,bool channel) { }; + Attribute(Attribute &a) { }; + ~Attribute() { }; + + void setUsed(bool set){ } + bool getUsed() { return true; } + void setIsChannel(bool set){ } + bool getIsChannel() { return false; } + + string getAsString(bool debug=false); + int getAsInt(); + bool getAsBool(); + double getAsFloat(); + ntlVec3d getAsVec3d(); + void getAsMat4Gfx(ntlMatrix4x4 *mat); + + AnimChannel getChannelInt(); + AnimChannel getChannelFloat(); + AnimChannel getChannelVec3d(); + AnimChannel getChannelSetVec3f(); + + string getCompleteString(); + void print(); + + protected: + + bool initChannel(int elemSize); +}; + + +// warning: DEPRECATED - replaced by elbeem API //! The list of configuration attributes class AttributeList { public: - //! Standard constructor - AttributeList(string name) : - mName(name), mAttrs() { }; - //! Destructor , delete all contained attribs + AttributeList(string name) { }; ~AttributeList(); - - /*! add an attribute to this list */ - void addAttr(string name, vector &value, int line, bool isChannel) { - if(exists(name)) delete mAttrs[name]; - mAttrs[name] = new Attribute(name,value,line, isChannel); - } - - /*! check if an attribute is set */ - bool exists(string name) { - if(mAttrs.find(name) == mAttrs.end()) return false; - return true; - } - - /*! get an attribute */ - Attribute *find(string name) { - if(mAttrs.find(name) == mAttrs.end()) { - errFatal("AttributeList::find","Invalid attribute '"< empty; - return new Attribute(name,empty, -1, 0); - } - return mAttrs[name]; - } - - //! set all params to used, for invisible objects + void addAttr(string name, vector &value, int line, bool isChannel) { } + bool exists(string name) { return false; } void setAllUsed(); - //! check if there were unknown params bool checkUnusedParams(); - - //! import attributes from other attribute list void import(AttributeList *oal); - - //! read attributes for object initialization int readInt(string name, int defaultValue, string source,string target, bool needed); bool readBool(string name, bool defaultValue, string source,string target, bool needed); double readFloat(string name, double defaultValue, string source,string target, bool needed); string readString(string name, string defaultValue, string source,string target, bool needed); ntlVec3d readVec3d(string name, ntlVec3d defaultValue, string source,string target, bool needed); void readMat4Gfx(string name, ntlMatrix4x4 defaultValue, string source,string target, bool needed, ntlMatrix4x4 *mat); - //! read attributes channels (attribute should be inited before) AnimChannel readChannelInt( string name, int defaultValue=0, string source=string("src"), string target=string("dst"), bool needed=false ); AnimChannel readChannelFloat( string name, double defaultValue=0, string source=string("src"), string target=string("dst"), bool needed=false ); AnimChannel readChannelVec3d( string name, ntlVec3d defaultValue=ntlVec3d(0.), string source=string("src"), string target=string("dst"), bool needed=false ); AnimChannel readChannelSetVec3f(string name, ntlSetVec3f defaultValue=ntlSetVec3f(0.), string source=string("src"), string target=string("dst"), bool needed=false ); - // channels with conversion AnimChannel readChannelVec3f( string name, ntlVec3f defaultValue=ntlVec3f(0.), string source=string("src"), string target=string("dst"), bool needed=false ); AnimChannel readChannelSinglePrecFloat( string name, float defaultValue=0., string source=string("src"), string target=string("dst"), bool needed=false ); - - //! set that a parameter can be given, and will be ignored... bool ignoreParameter(string name, string source); - - //! debug function, prints all attribs void print(); - protected: - - /*! attribute name (form config file) */ - string mName; - - /*! the global attribute storage */ - map mAttrs; - }; ntlVec3f channelFindMaxVf (AnimChannel channel); @@ -297,6 +203,14 @@ bool channelSimplifyi (AnimChannel &channel); bool channelSimplifyf (AnimChannel &channel); bool channelSimplifyd (AnimChannel &channel); +//! output channel values? on=1/off=0 +#define DEBUG_PCHANNELS 0 + +//! debug function, prints to stdout if DEBUG_PCHANNELS flag is enabled, used in constructors +template +void AnimChannel::debugPrintChannel() { } + + #define NTL_ATTRIBUTES_H #endif diff --git a/intern/elbeem/intern/elbeem.cpp b/intern/elbeem/intern/elbeem.cpp index d4242da4d5e..b2779f51c3b 100644 --- a/intern/elbeem/intern/elbeem.cpp +++ b/intern/elbeem/intern/elbeem.cpp @@ -3,7 +3,7 @@ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method * All code distributed as part of El'Beem is covered by the version 2 of the * GNU General Public License. See the file COPYING for details. - * Copyright 2003-2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Main program functions */ @@ -64,10 +64,11 @@ void elbeemResetSettings(elbeemSimulationSettings *set) { set->channelSizeGravity=0; set->channelGravity=NULL; - set->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP; - set->obstaclePartslip= 0.; + set->domainobsType= FLUIDSIM_OBSTACLE_NOSLIP; + set->domainobsPartslip= 0.; set->generateVertexVectors = 0; set->surfaceSmoothing = 1.; + set->surfaceSubdivs = 1; set->farFieldSize = 0.; set->runsimCallback = NULL; @@ -83,6 +84,7 @@ extern "C" int elbeemInit() { setElbeemState( SIMWORLD_INITIALIZING ); setElbeemErrorString("[none]"); + resetGlobalColorSetting(); elbeemCheckDebugEnv(); debMsgStd("performElbeemSimulation",DM_NOTIFY,"El'Beem Simulation Init Start as Plugin, debugLevel:"<obstacleType= FLUIDSIM_OBSTACLE_NOSLIP; mesh->obstaclePartslip= 0.; + mesh->obstacleImpactFactor= 1.; mesh->volumeInitType= OB_VOLUMEINIT_VOLUME; @@ -193,6 +196,7 @@ int elbeemAddMesh(elbeemMesh *mesh) { obj->setGeoInitIntersect(true); obj->setGeoInitType(initType); obj->setGeoPartSlipValue(mesh->obstaclePartslip); + obj->setGeoImpactFactor(mesh->obstacleImpactFactor); if((mesh->volumeInitTypevolumeInitType>VOLUMEINIT_BOTH)) mesh->volumeInitType = VOLUMEINIT_VOLUME; obj->setVolumeInit(mesh->volumeInitType); // use channel instead, obj->setInitialVelocity( ntlVec3Gfx(mesh->iniVelocity[0], mesh->iniVelocity[1], mesh->iniVelocity[2]) ); diff --git a/intern/elbeem/intern/elbeem.h b/intern/elbeem/intern/elbeem.h index ea21e807a6f..b3feda8bbe8 100644 --- a/intern/elbeem/intern/elbeem.h +++ b/intern/elbeem/intern/elbeem.h @@ -3,7 +3,7 @@ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method * All code distributed as part of El'Beem is covered by the version 2 of the * GNU General Public License. See the file COPYING for details. - * Copyright 2003-2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * API header */ @@ -72,13 +72,14 @@ typedef struct elbeemSimulationSettings { float *channelGravity; // vector /* boundary types and settings for domain walls */ - short obstacleType; - float obstaclePartslip; + short domainobsType; + float domainobsPartslip; /* generate speed vectors for vertices (e.g. for image based motion blur)*/ short generateVertexVectors; /* strength of surface smoothing */ float surfaceSmoothing; - // TODO add surf gen flags + /* no. of surface subdivisions */ + int surfaceSubdivs; /* global transformation to apply to fluidsim mesh */ float surfaceTrafo[4*4]; @@ -147,6 +148,8 @@ typedef struct elbeemMesh { /* boundary types and settings */ short obstacleType; float obstaclePartslip; + /* amount of force transfer from fluid to obj, 0=off, 1=normal */ + float obstacleImpactFactor; /* init volume, shell or both? use OB_VOLUMEINIT_xxx defines above */ short volumeInitType; @@ -230,5 +233,8 @@ double elbeemEstimateMemreq(int res, #define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6)) #define FGI_MBNDOUTFLOW (1<<(FGI_FLAGSTART+ 7)) +// all boundary types at once +#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW ) + #endif // ELBEEM_API_H diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp index 283fd17eb94..fdef04dc5c6 100644 --- a/intern/elbeem/intern/isosurface.cpp +++ b/intern/elbeem/intern/isosurface.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Marching Cubes surface mesh generation * @@ -9,12 +9,10 @@ #include "isosurface.h" #include "mcubes_tables.h" -#include "ntl_ray.h" - +#include "particletracer.h" #include #include - // sirdude fix for solaris #if !defined(linux) && (defined (__sparc) || defined (__sparc__)) #include @@ -30,14 +28,17 @@ IsoSurface::IsoSurface(double iso) : mpData(NULL), mIsoValue( iso ), mPoints(), + mUseFullEdgeArrays(false), mpEdgeVerticesX(NULL), mpEdgeVerticesY(NULL), mpEdgeVerticesZ(NULL), + mEdgeArSize(-1), mIndices(), mStart(0.0), mEnd(0.0), mDomainExtent(0.0), mInitDone(false), mSmoothSurface(0.0), mSmoothNormals(0.0), mAcrossEdge(), mAdjacentFaces(), - mCutoff(-1), mCutArray(NULL),// off by default + mCutoff(-1), mCutArray(NULL), // off by default + mpIsoParts(NULL), mPartSize(0.), mSubdivs(0), mFlagCnt(1), mSCrad1(0.), mSCrad2(0.), mSCcenter(0.) { @@ -49,6 +50,10 @@ IsoSurface::IsoSurface(double iso) : *****************************************************************************/ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent) { + // range 1-10 (max due to subd array in triangulate) + if(mSubdivs<1) mSubdivs=1; + if(mSubdivs>10) mSubdivs=10; + // init solver and size mSizex = setx; mSizey = sety; @@ -74,14 +79,27 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e for(int i=0;i0) continue; // */ + value[0] = *getData(i ,j ,k ); + value[1] = *getData(i+1,j ,k ); + value[2] = *getData(i+1,j+1,k ); + value[3] = *getData(i ,j+1,k ); + value[4] = *getData(i ,j ,k+1); + value[5] = *getData(i+1,j ,k+1); + value[6] = *getData(i+1,j+1,k+1); + value[7] = *getData(i ,j+1,k+1); - // check intersections of isosurface with edges, and calculate cubie index - cubeIndex = 0; - if (value[0] < mIsoValue) cubeIndex |= 1; - if (value[1] < mIsoValue) cubeIndex |= 2; - if (value[2] < mIsoValue) cubeIndex |= 4; - if (value[3] < mIsoValue) cubeIndex |= 8; - if (value[4] < mIsoValue) cubeIndex |= 16; - if (value[5] < mIsoValue) cubeIndex |= 32; - if (value[6] < mIsoValue) cubeIndex |= 64; - if (value[7] < mIsoValue) cubeIndex |= 128; + // check intersections of isosurface with edges, and calculate cubie index + cubeIndex = 0; + if (value[0] < mIsoValue) cubeIndex |= 1; + if (value[1] < mIsoValue) cubeIndex |= 2; + if (value[2] < mIsoValue) cubeIndex |= 4; + if (value[3] < mIsoValue) cubeIndex |= 8; + if (value[4] < mIsoValue) cubeIndex |= 16; + if (value[5] < mIsoValue) cubeIndex |= 32; + if (value[6] < mIsoValue) cubeIndex |= 64; + if (value[7] < mIsoValue) cubeIndex |= 128; - // No triangles to generate? - if (mcEdgeTable[cubeIndex] == 0) { - continue; - } - - // where to look up if this point already exists - eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; - eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; - eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+0) ]; - eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; - - eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, k+1) ]; - eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, k+1) ]; - eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, k+1) ]; - eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, k+1) ]; - - eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; - eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, k+0) ]; - eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, k+0) ]; - eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, k+0) ]; - - // grid positions - pos[0] = ntlVec3Gfx(px ,py ,pz); - pos[1] = ntlVec3Gfx(px+gsx,py ,pz); - pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); - pos[3] = ntlVec3Gfx(px ,py+gsy,pz); - pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); - pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); - pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); - pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); - - // check all edges - for(int e=0;e<12;e++) { - - if (mcEdgeTable[cubeIndex] & (1<0) && (kmSizex-2-coAdd-mCutoff) || - (j>mSizey-2-coAdd-mCutoff) ) { - if(mCutArray) { - if(k < mCutArray[j*this->mSizex+i]) continue; - } else { + // No triangles to generate? + if (mcEdgeTable[cubeIndex] == 0) { continue; } - } - // Create the triangles... - for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) { - //errMsg("MC","tri "<0) && (kmSizex-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) { + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+0] ] ); + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+1] ] ); + mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] ); + } + + }//i + }// j + + // copy edge arrays + if(!mUseFullEdgeArrays) { + for(int j=0;j<(mSizey-0);j++) + for(int i=0;i<(mSizex-0);i++) { + //int edgek = 0; + const int dst = ISOLEVEL_INDEX( i+0, j+0, 0); + const int src = ISOLEVEL_INDEX( i+0, j+0, 1); + mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ]; + mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; + mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ]; + mpEdgeVerticesX[ src ]=-1; + mpEdgeVerticesY[ src ]=-1; + mpEdgeVerticesZ[ src ]=-1; + } + } // */ + + } // k + + // precalculate normals using an approximation of the scalar field gradient + for(int ni=0;ni<(int)mPoints.size();ni++) { normalize( mPoints[ni].n ); } + + } else { // subdivs + +#define EDGEAR_INDEX(Ai,Aj,Ak, Bi,Bj) ((mSizex*mSizey*mSubdivs*mSubdivs*(Ak))+\ + (mSizex*mSubdivs*((Aj)*mSubdivs+(Bj)))+((Ai)*mSubdivs)+(Bi)) + +#define ISOTRILININT(fi,fj,fk) ( \ + (1.-(fi))*(1.-(fj))*(1.-(fk))*orgval[0] + \ + ( (fi))*(1.-(fj))*(1.-(fk))*orgval[1] + \ + ( (fi))*( (fj))*(1.-(fk))*orgval[2] + \ + (1.-(fi))*( (fj))*(1.-(fk))*orgval[3] + \ + (1.-(fi))*(1.-(fj))*( (fk))*orgval[4] + \ + ( (fi))*(1.-(fj))*( (fk))*orgval[5] + \ + ( (fi))*( (fj))*( (fk))*orgval[6] + \ + (1.-(fi))*( (fj))*( (fk))*orgval[7] ) + + // use subdivisions + gfxReal subdfac = 1./(gfxReal)(mSubdivs); + gfxReal orgGsx = gsx; + gfxReal orgGsy = gsy; + gfxReal orgGsz = gsz; + gsx *= subdfac; + gsy *= subdfac; + gsz *= subdfac; + if(mUseFullEdgeArrays) { + errMsg("IsoSurface::triangulate","Disabling mUseFullEdgeArrays!"); + } + + // subdiv local arrays + gfxReal orgval[8]; + gfxReal subdAr[2][11][11]; // max 10 subdivs! + ParticleObject* *arppnt = new ParticleObject*[mSizez*mSizey*mSizex]; + + // construct pointers + // part test + int pInUse = 0; + int pUsedTest = 0; + // reset particles + // reset list array + for(int k=0;k<(mSizez);k++) + for(int j=0;j<(mSizey);j++) + for(int i=0;i<(mSizex);i++) { + arppnt[ISOLEVEL_INDEX(i,j,k)] = NULL; + } + if(mpIsoParts) { + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { + if( (*pit).getActive()==false ) continue; + if( (*pit).getType()!=PART_DROP) continue; + (*pit).setNext(NULL); + } + // build per node lists + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { + if( (*pit).getActive()==false ) continue; + if( (*pit).getType()!=PART_DROP) continue; + // check lifetime ignored here + ParticleObject *p = &(*pit); + const ntlVec3Gfx ppos = p->getPos(); + const int pi= (int)round(ppos[0])+0; + const int pj= (int)round(ppos[1])+0; + int pk= (int)round(ppos[2])+0;// no offset necessary + // 2d should be handled by solver. if(LBMDIM==2) { pk = 0; } + + if(pi<0) continue; + if(pj<0) continue; + if(pk<0) continue; + if(pi>mSizex-1) continue; + if(pj>mSizey-1) continue; + if(pk>mSizez-1) continue; + ParticleObject* &pnt = arppnt[ISOLEVEL_INDEX(pi,pj,pk)]; + if(pnt) { + // append + ParticleObject* listpnt = pnt; + while(listpnt) { + if(!listpnt->getNext()) { + listpnt->setNext(p); listpnt = NULL; + } else { + listpnt = listpnt->getNext(); + } + } + } else { + // start new list + pnt = p; + } + pInUse++; + } + } // mpIsoParts + + debMsgStd("IsoSurface::triangulate",DM_MSG,"Starting. Parts in use:"<=mSizez-1) continue; + for(int poj=-poDistOffset; poj<1+poDistOffset; poj++) { + if(j+poj<0) continue; + if(j+poj>=mSizey-1) continue; + for(int poi=-poDistOffset; poi<1+poDistOffset; poi++) { + if(i+poi<0) continue; + if(i+poi>=mSizex-1) continue; + ParticleObject *p; + p = arppnt[ISOLEVEL_INDEX(i+poi,j+poj,k+pok)]; + while(p) { // */ + /* + for(vector::iterator pit= mpIsoParts->getParticlesBegin(); + pit!= mpIsoParts->getParticlesEnd(); pit++) { { { { + // debug test! , full list slow! + if(( (*pit).getActive()==false ) || ( (*pit).getType()!=PART_DROP)) continue; + ParticleObject *p; + p = &(*pit); // */ + + pUsedTest++; + ntlVec3Gfx ppos = p->getPos(); + const int spi= (int)round( (ppos[0]+1.-(gfxReal)i) *(gfxReal)mSubdivs-1.5); + const int spj= (int)round( (ppos[1]+1.-(gfxReal)j) *(gfxReal)mSubdivs-1.5); + const int spk= (int)round( (ppos[2]+1.-(gfxReal)k) *(gfxReal)mSubdivs-1.5)-sdkOffset; // why -2? + // 2d should be handled by solver. if(LBMDIM==2) { spk = 0; } + + gfxReal pfLen = p->getSize()*1.5*mPartSize; // test, was 1.1 + const gfxReal minPfLen = subdfac*0.8; + if(pfLengetSize()<<" ps"< 1) { continue; } // */ + for(int swj=-icellpsize; swj<=icellpsize; swj++) { + if(spj+swj< 0) { continue; } + if(spj+swj>mSubdivs+0) { continue; } // */ + for(int swi=-icellpsize; swi<=icellpsize; swi++) { + if(spi+swi< 0) { continue; } + if(spi+swi>mSubdivs+0) { continue; } // */ + ntlVec3Gfx cellp = ntlVec3Gfx( + (1.5+(gfxReal)(spi+swi)) *subdfac + (gfxReal)(i-1), + (1.5+(gfxReal)(spj+swj)) *subdfac + (gfxReal)(j-1), + (1.5+(gfxReal)(spk+swk)+sdkOffset) *subdfac + (gfxReal)(k-1) + ); + //if(swi==0 && swj==0 && swk==0) subdAr[spk][spj][spi] = 1.; // DEBUG + // clip domain boundaries again + if(cellp[0]<1.) { continue; } + if(cellp[1]<1.) { continue; } + if(cellp[2]<1.) { continue; } + if(cellp[0]>(gfxReal)mSizex-3.) { continue; } + if(cellp[1]>(gfxReal)mSizey-3.) { continue; } + if(cellp[2]>(gfxReal)mSizez-3.) { continue; } + gfxReal len = norm(cellp-ppos); + gfxReal isoadd = 0.; + const gfxReal baseIsoVal = mIsoValue*1.1; + if(len1.) { continue; } + subdAr[spk+swk][spj+swj][spi+swi] = arval + isoadd; + } } } + + p = p->getNext(); + } + } } } // poDist loops */ + + py = mStart[1]+(((double)j-0.5)*orgGsy)-gsy; + for(int sj=0;sj 0) { + + // where to look up if this point already exists + const int edgek = 0; + const int baseIn = EDGEAR_INDEX( i+0, j+0, edgek+0, si,sj); + eVert[ 0] = &mpEdgeVerticesX[ baseIn ]; + eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ]; + eVert[ 2] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; + eVert[ 3] = &mpEdgeVerticesY[ baseIn ]; + + eVert[ 4] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; + eVert[ 5] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+1,sj+0) ]; // with subdivs + eVert[ 6] = &mpEdgeVerticesX[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+1) ]; + eVert[ 7] = &mpEdgeVerticesY[ EDGEAR_INDEX( i, j, edgek+1, si+0,sj+0) ]; + + eVert[ 8] = &mpEdgeVerticesZ[ baseIn ]; + eVert[ 9] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+0) ]; // with subdivs + eVert[10] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+1,sj+1) ]; + eVert[11] = &mpEdgeVerticesZ[ EDGEAR_INDEX( i, j, edgek+0, si+0,sj+1) ]; + + // grid positions + pos[0] = ntlVec3Gfx(px ,py ,pz); + pos[1] = ntlVec3Gfx(px+gsx,py ,pz); + pos[2] = ntlVec3Gfx(px+gsx,py+gsy,pz); // with subdivs + pos[3] = ntlVec3Gfx(px ,py+gsy,pz); + pos[4] = ntlVec3Gfx(px ,py ,pz+gsz); + pos[5] = ntlVec3Gfx(px+gsx,py ,pz+gsz); + pos[6] = ntlVec3Gfx(px+gsx,py+gsy,pz+gsz); // with subdivs + pos[7] = ntlVec3Gfx(px ,py+gsy,pz+gsz); + + // check all edges + for(int e=0;e<12;e++) { + if (mcEdgeTable[cubeIndex] & (1<0) { + //smoSubdfac = 1./(float)(mSubdivs); + smoSubdfac = powf(0.55,(float)mSubdivs); // slightly stronger + } + if(mSmoothSurface>0. || mSmoothNormals>0.) debMsgStd("IsoSurface::triangulate",DM_MSG,"Smoothing...",10); + if(mSmoothSurface>0.0) { + smoothSurface(mSmoothSurface*smoSubdfac, (mSmoothNormals<=0.0) ); + } + if(mSmoothNormals>0.0) { + smoothNormals(mSmoothNormals*smoSubdfac); } - if(mSmoothSurface>0.0) { smoothSurface(mSmoothSurface, (mSmoothNormals<=0.0) ); } - if(mSmoothNormals>0.0) { smoothNormals(mSmoothNormals); } myTime_t tritimeend = getTime(); - debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<getNumParticles(), 10); } @@ -405,12 +750,8 @@ inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) { /****************************************************************************** * - * Surface improvement - * makes use of trimesh2 library - * http://www.cs.princeton.edu/gfx/proj/trimesh2/ - * - * Copyright (c) 2004 Szymon Rusinkiewicz. - * see COPYING_trimesh2 + * Surface improvement, inspired by trimesh2 library + * (http://www.cs.princeton.edu/gfx/proj/trimesh2/) * *****************************************************************************/ @@ -420,9 +761,39 @@ void IsoSurface::setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc) { mSCcenter = mscc; } +// compute normals for all generated triangles +void IsoSurface::computeNormals() { + for(int i=0;i<(int)mPoints.size();i++) { + mPoints[i].n = ntlVec3Gfx(0.); + } + + for(int i=0;i<(int)mIndices.size();i+=3) { + const int t1 = mIndices[i]; + const int t2 = mIndices[i+1]; + const int t3 = mIndices[i+2]; + const ntlVec3Gfx p1 = mPoints[t1].v; + const ntlVec3Gfx p2 = mPoints[t2].v; + const ntlVec3Gfx p3 = mPoints[t3].v; + const ntlVec3Gfx n1=p1-p2; + const ntlVec3Gfx n2=p2-p3; + const ntlVec3Gfx n3=p3-p1; + const gfxReal len1 = normNoSqrt(n1); + const gfxReal len2 = normNoSqrt(n2); + const gfxReal len3 = normNoSqrt(n3); + const ntlVec3Gfx norm = cross(n1,n2); + mPoints[t1].n += norm * (1./(len1*len3)); + mPoints[t2].n += norm * (1./(len1*len2)); + mPoints[t3].n += norm * (1./(len2*len3)); + } + + for(int i=0;i<(int)mPoints.size();i++) { + normalize(mPoints[i].n); + } +} + // Diffuse a vector field at 1 vertex, weighted by // a gaussian of width 1/sqrt(invsigma2) -bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &flt) +bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int src, float invsigma2, ntlVec3Gfx &target) { if((neighbors[src].size()<1) || (pointareas[src]<=0.0)) return 0; const ntlVec3Gfx srcp = mPoints[src].v; @@ -431,31 +802,23 @@ bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, i ntlVec3Gfx dp = mSCcenter-srcp; dp[2] = 0.0; // only xy-plane float rd = normNoSqrt(dp); if(rd > mSCrad2) { - //errMsg("TRi","p"< mSCrad1) { // optimize? float org = 1.0/sqrt(invsigma2); org *= (1.0- (rd-mSCrad1) / (mSCrad2-mSCrad1)); invsigma2 = 1.0/(org*org); - //flt = ntlVec3Gfx((rd-mSCrad1) / (mSCrad2-mSCrad1)); return 1; //errMsg("TRi","p"< mDboundary = neighbors[src]; mDboundary = neighbors[src]; while (!mDboundary.empty()) { const int bbn = mDboundary.back(); @@ -465,33 +828,17 @@ bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, i // normal check const float nvdot = dot(srcn, mPoints[bbn].n); // faster than before d2 calc? - if(nvdot <= 0.0f) continue; // faster than before d2 calc? + if(nvdot <= 0.0f) continue; // gaussian weight of width 1/sqrt(invsigma2) const float d2 = invsigma2 * normNoSqrt(mPoints[bbn].v - srcp); - if(d2 >= 9.0f) continue; // 25 also possible , slower - //if(dot(srcn, mPoints[bbn].n) <= 0.0f) continue; // faster than before d2 calc? + if(d2 >= 9.0f) continue; - //float w = (d2 >= 9.0f) ? 0.0f : exp(-0.5f*d2); - //float w = expf(-0.5f*d2); -#if 0 - float w=1.0; - // Downweight things pointing in different directions - w *= nvdot; //dot(srcn , mPoints[bbn].n); - // Surface area "belonging" to each point - w *= pointareas[bbn]; + // aggressive smoothing factor + float smstr = nvdot * pointareas[bbn]; // Accumulate weight times field at neighbor - flt += *(field+pointerScale*bbn)*w; - sum_w += w; - // */ -#else - // more aggressive smoothing with: float w=1.0; - float w=nvdot * pointareas[bbn]; - // Accumulate weight times field at neighbor - flt += *(field+pointerScale*bbn)*w; - sum_w += w; -#endif - // */ + target += *(field+pointerScale*bbn)*smstr; + smstrSum += smstr; for(int i = 0; i < (int)neighbors[bbn].size(); i++) { const int nn = neighbors[bbn][i]; @@ -499,17 +846,12 @@ bool IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, i mDboundary.push_back(nn); } } - flt /= sum_w; + target /= smstrSum; return 1; } -// REF -// TestData::getTriangles message: Time for surface generation:3.75s, S(0.0390625,0.1171875) - // ntlWorld::ntlWorld message: Time for start-sims:0s - // TestData::getTriangles message: Time for surface generation:3.69s, S(0.0390625,0.1171875) - - +// perform smoothing of the surface (and possible normals) void IsoSurface::smoothSurface(float sigma, bool normSmooth) { int nv = mPoints.size(); @@ -670,8 +1012,8 @@ void IsoSurface::smoothSurface(float sigma, bool normSmooth) //errMsg("SMSURF","done v:"<10) s=10; + mSubdivs = s; } + int getSubdivs() { return mSubdivs;} + /*! set full edge settings, this has to be done before init! */ + void setUseFulledgeArrays(bool set) { + if(mInitDone) errFatal("IsoSurface::setUseFulledgeArrays","Changing usefulledge after init!", SIMWORLD_INITERROR); + mUseFullEdgeArrays = set;} + protected: /* variables ... */ @@ -69,10 +84,13 @@ class IsoSurface : //! Store all the triangles vertices vector mPoints; + //! use full arrays? (not for farfield) + bool mUseFullEdgeArrays; //! Store indices of calculated points along the cubie edges int *mpEdgeVerticesX; int *mpEdgeVerticesY; int *mpEdgeVerticesZ; + int mEdgeArSize; //! vector for all the triangles (stored as 3 indices) @@ -100,6 +118,12 @@ class IsoSurface : int mCutoff; //! cutoff heigh values int *mCutArray; + //! particle pointer + ParticleTracer *mpIsoParts; + //! particle size + float mPartSize; + //! no of subdivisions + int mSubdivs; //! trimesh vars vector flags; @@ -186,6 +210,7 @@ class IsoSurface : void setSmoothRad(float radi1, float radi2, ntlVec3Gfx mscc); void smoothSurface(float val, bool smoothNorm); void smoothNormals(float val); + void computeNormals(); protected: diff --git a/intern/elbeem/intern/loop_tools.h b/intern/elbeem/intern/loop_tools.h index 927263b1204..3c15a609210 100644 --- a/intern/elbeem/intern/loop_tools.h +++ b/intern/elbeem/intern/loop_tools.h @@ -30,6 +30,7 @@ #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz); #define LIST_EMPTY(x) mListEmpty.push_back( x ); #define LIST_FULL(x) mListFull.push_back( x ); +#define FSGR_ADDPART(x) mpParticles->addFullParticle( x ); // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #define GRID_REGION_START() \ @@ -90,6 +91,9 @@ int i=0; \ ADVANCE_POINTERS(2*gridLoopBound); \ } /* j */ \ + /* COMPRESSGRIDS!=1 */ \ + /* int i=0; */ \ + /* ADVANCE_POINTERS(mLevel[lev].lSizex*2); */ \ } /* all cell loop k,j,i */ \ if(doReduce) { } /* dummy remove warning */ \ } /* main_region */ \ diff --git a/intern/elbeem/intern/ntl_blenderdumper.cpp b/intern/elbeem/intern/ntl_blenderdumper.cpp index 1bdea59f210..b1fece25890 100644 --- a/intern/elbeem/intern/ntl_blenderdumper.cpp +++ b/intern/elbeem/intern/ntl_blenderdumper.cpp @@ -1,13 +1,12 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Replaces std. raytracer, and only dumps time dep. objects to disc * *****************************************************************************/ - #include #include @@ -44,6 +43,10 @@ ntlBlenderDumper::~ntlBlenderDumper() debMsgStd("ntlBlenderDumper",DM_NOTIFY, "ntlBlenderDumper done", 10); } +// required globals +extern bool glob_mpactive; +extern int glob_mpnum, glob_mpindex; + /****************************************************************************** * Only dump time dep. objects to file *****************************************************************************/ @@ -52,7 +55,8 @@ int ntlBlenderDumper::renderScene( void ) char nrStr[5]; /* nr conversion */ ntlRenderGlobals *glob = mpGlob; ntlScene *scene = mpGlob->getSimScene(); - bool debugOut = true; + bool debugOut = false; + bool debugRender = false; #if ELBEEM_PLUGIN==1 debugOut = false; #endif // ELBEEM_PLUGIN==1 @@ -63,8 +67,6 @@ int ntlBlenderDumper::renderScene( void ) if(debugOut) debMsgStd("ntlBlenderDumper::renderScene",DM_NOTIFY,"Dumping geometry data", 1); long startTime = getTime(); - - /* check if picture already exists... */ snprintf(nrStr, 5, "%04d", glob->getAniCount() ); // local scene vars @@ -72,7 +74,7 @@ int ntlBlenderDumper::renderScene( void ) vector Vertices; vector VertNormals; - /* init geometry array, first all standard objects */ + // check geo objects int idCnt = 0; // give IDs to objects for (vector::iterator iter = scene->getGeoClasses()->begin(); iter != scene->getGeoClasses()->end(); iter++) { @@ -80,8 +82,7 @@ int ntlBlenderDumper::renderScene( void ) int tid = (*iter)->getTypeId(); if(tid & GEOCLASSTID_OBJECT) { - // normal geom. objects dont change... -> ignore - //if(buildInfo) debMsgStd("ntlBlenderDumper::BuildScene",DM_MSG,"added GeoObj "<getName(), 8 ); + // normal geom. objects -> ignore } if(tid & GEOCLASSTID_SHADER) { ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter); //dynamic_cast(*iter); @@ -105,6 +106,11 @@ int ntlBlenderDumper::renderScene( void ) isPreview = true; } if(!doDump) continue; + + // dont quit, some objects need notifyOfDump call + if((glob_mpactive) && (glob_mpindex>0)) { + continue; //return 0; + } // only dump geo shader objects Triangles.clear(); @@ -112,8 +118,12 @@ int ntlBlenderDumper::renderScene( void ) VertNormals.clear(); (*siter)->initialize( mpGlob ); (*siter)->getTriangles(this->mSimulationTime, &Triangles, &Vertices, &VertNormals, idCnt); - idCnt ++; + + // WARNING - this is dirty, but simobjs are the only geoshaders right now + SimulationObject *sim = (SimulationObject *)geoshad; + LbmSolverInterface *lbm = sim->getSolver(); + // always dump mesh, even empty ones... @@ -127,9 +137,6 @@ int ntlBlenderDumper::renderScene( void ) gzFile gzf; // output velocities if desired - // WARNING - this is dirty, but simobjs are the only geoshaders right now - SimulationObject *sim = (SimulationObject *)geoshad; - LbmSolverInterface *lbm = sim->getSolver(); if((!isPreview) && (lbm->getDumpVelocities())) { std::ostringstream bvelfilename; bvelfilename << boutfilename.str(); @@ -155,24 +162,38 @@ int ntlBlenderDumper::renderScene( void ) // compress all bobj's boutfilename << ".bobj.gz"; - //if(isPreview) { } else { } gzf = gzopen(boutfilename.str().c_str(), "wb1"); // wb9 is slow for large meshes! if (!gzf) { errMsg("ntlBlenderDumper::renderScene","Unable to open output '"< *trafo; trafo = lbm->getDomainTrafo(); if(trafo) { - //trafo->initArrayCheck(ettings->surfaceTrafo); - //errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) ); // transform into source space for(size_t i=0; i rottrafo; + rottrafo.initId(); + if(lbm->getDomainTrafo()) { + // dont modifiy original! + rottrafo = *lbm->getDomainTrafo(); + ntlVec3Gfx rTrans,rScale,rRot,rShear; + rottrafo.decompose(rTrans,rScale,rRot,rShear); + rottrafo.initRotationXYZ(rRot[0],rRot[1],rRot[2]); + // only rotate here... + for(size_t i=0; i0) { if(debugOut) debMsgStd("ntlBlenderDumper::renderScene",DM_MSG,"Objects dumped: "<0)) { + // ok, nothing to do anyway... + } else { + errFatal("ntlBlenderDumper::renderScene","No objects to dump! Aborting...",SIMWORLD_INITERROR); + return 1; + } } - /* next frame */ - //glob->setAniCount( glob->getAniCount() +1 ); + // debug timing long stopTime = getTime(); debMsgStd("ntlBlenderDumper::renderScene",DM_MSG,"Scene #"<1) || (mcScale.accessValues().size()>1) || (mcGeoActive.accessValues().size()>1) + // mcGeoImpactFactor only needed when moving || (mcInitialVelocity.accessValues().size()>1) ) { mIsAnimated = true; @@ -71,6 +75,7 @@ bool ntlGeometryObject::checkIsAnimated() { if(mGeoInitType==FGI_FLUID) { mIsAnimated=false; } + //errMsg("ntlGeometryObject::checkIsAnimated","obj="<readInt("geoinit_volumeinit", mVolumeInit,"ntlGeometryObject", "mVolumeInit", false); if((mVolumeInitVOLUMEINIT_BOTH)) mVolumeInit = VOLUMEINIT_VOLUME; + // moving obs correction factor + float impactfactor=1.; + impactfactor = (float)mpAttrs->readFloat("impactfactor", impactfactor,"ntlGeometryObject", "impactfactor", false); + if(getAttributeList()->exists("impactfactor") || (!mcGeoImpactFactor.isInited()) ) { + mcGeoImpactFactor = mpAttrs->readChannelSinglePrecFloat("impactfactor"); + } + // override cfg types mVisible = mpAttrs->readBool("visible", mVisible,"ntlGeometryObject", "mVisible", false); mReceiveShadows = mpAttrs->readBool("recv_shad", mReceiveShadows,"ntlGeometryObject", "mReceiveShadows", false); @@ -162,12 +174,12 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob) } float geoactive=1.; - geoactive = mpAttrs->readFloat("geoactive", geoactive,"ntlGeometryObject", "geoactive", false); + geoactive = (float)mpAttrs->readFloat("geoactive", geoactive,"ntlGeometryObject", "geoactive", false); if(getAttributeList()->exists("geoactive") || (!mcGeoActive.isInited()) ) { - mcGeoActive = mpAttrs->readChannelFloat("geoactive"); + mcGeoActive = mpAttrs->readChannelSinglePrecFloat("geoactive"); } // always use channel - if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel(geoactive); } + if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel(geoactive); } checkIsAnimated(); @@ -309,12 +321,12 @@ void ntlGeometryObject::sceneAddTriangleNoVert(int *trips, (dst) = AnimChannel< ntlVec3Gfx >(vals,time); #define ADD_CHANNEL_FLOAT(dst,nvals,val) \ - valsd.clear(); time.clear(); elbeemSimplifyChannelFloat(val,&nvals); \ + valsfloat.clear(); time.clear(); elbeemSimplifyChannelFloat(val,&nvals); \ for(int i=0; i<(nvals); i++) { \ - valsd.push_back( (val)[i*2+0] ); \ + valsfloat.push_back( (val)[i*2+0] ); \ time.push_back( (val)[i*2+1] ); \ } \ - (dst) = AnimChannel< double >(valsd,time); + (dst) = AnimChannel< float >(valsfloat,time); void ntlGeometryObject::initChannels( int nTrans, float *trans, int nRot, float *rot, int nScale, float *scale, @@ -324,7 +336,7 @@ void ntlGeometryObject::initChannels( if(debugInitc) { debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"nt:"< vals; - vector valsd; + vector valsfloat; vector time; if((trans)&&(nTrans>0)) { ADD_CHANNEL_VEC(mcTrans, nTrans, trans); } if((rot)&&(nRot>0)) { ADD_CHANNEL_VEC(mcRot, nRot, rot); } @@ -383,7 +395,7 @@ void ntlGeometryObject::applyTransformation(double t, vector *verts, ntlMat4Gfx rotMat; rotMat.initRotationXYZ(rot[0],rot[1],rot[2]); pos += mInitialPos; - //errMsg("ntlGeometryObject::applyTransformation","obj="< *verts, } } else { // not animated, cached points were already returned - //errMsg ("ntlGeometryObject::applyTransformation","Object "< 0) { for(int u=0; u<=divs1; u++) { for(int v=0; v<=divs2; v++) { - const gfxReal uf = (gfxReal)(u+0.25) / (gfxReal)(divs1+0.0); - const gfxReal vf = (gfxReal)(v+0.25) / (gfxReal)(divs2+0.0); + const gfxReal uf = (gfxReal)(u+TRI_UVOFFSET) / (gfxReal)(divs1+0.0); + const gfxReal vf = (gfxReal)(v+TRI_UVOFFSET) / (gfxReal)(divs2+0.0); if(uf+vf>1.0) continue; countp+=2; } @@ -500,7 +512,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) { } mMovPoints.push_back(p); mMovNormals.push_back(n); - //errMsg("ntlGeometryObject::initMovingPoints","std"< 0) { for(int u=0; u<=divs1; u++) { for(int v=0; v<=divs2; v++) { - const gfxReal uf = (gfxReal)(u+0.25) / (gfxReal)(divs1+0.0); - const gfxReal vf = (gfxReal)(v+0.25) / (gfxReal)(divs2+0.0); + const gfxReal uf = (gfxReal)(u+TRI_UVOFFSET) / (gfxReal)(divs1+0.0); + const gfxReal vf = (gfxReal)(v+TRI_UVOFFSET) / (gfxReal)(divs2+0.0); if(uf+vf>1.0) continue; ntlVec3Gfx p = vertices[ trips[0] ] * (1.0-uf-vf)+ @@ -653,8 +665,8 @@ void ntlGeometryObject::initMovingPointsAnim( //errMsg("ntlGeometryObject::initMovingPointsAnim","Tri1 "<1.0) continue; ntlVec3Gfx srcp = srcvertices[ srctrips[0] ] * (1.0-uf-vf)+ @@ -722,7 +734,7 @@ ntlVec3Gfx ntlGeometryObject::calculateMaxVel(double t1, double t2) { applyTransformation(t2,&verts2,NULL, 0,verts2.size(), true); vel = (verts2[0]-verts1[0]); // /(t2-t1); - errMsg("ntlGeometryObject::calculateMaxVel","t1="<size()<<" t"<size();i++) { SimulationObject *sim = (*mpSims)[i]; @@ -496,6 +500,8 @@ void ntlWorld::singleStepSims(double targetTime) { +extern bool glob_mpactive; +extern int glob_mpindex; /****************************************************************************** * Render the current scene @@ -511,11 +517,19 @@ int ntlWorld::renderScene( void ) myTime_t rendStart,rendEnd; // measure user rendering time glob = mpGlob; + // deactivate for all with index!=0 + if((glob_mpactive)&&(glob_mpindex>0)) return(0); + /* check if picture already exists... */ if(!glob->getSingleFrameMode() ) { snprintf(nrStr, 5, "%04d", glob->getAniCount() ); - //outfilename << glob->getOutFilename() <<"_" << nrStr << ".ppm"; - outfn_conv << glob->getOutFilename() <<"_" << nrStr << ".png"; + + if(glob_mpactive) { + outfn_conv << glob->getOutFilename() <<"_"<getOutFilename() <<"_" << nrStr << ".png"; + } //if((mpGlob->getDisplayMode() == DM_RAY)&&(mpGlob->getFrameSkip())) { if(mpGlob->getFrameSkip()) { @@ -823,11 +837,7 @@ int ntlWorld::renderScene( void ) } for(int i = 0; i < h; i++) rows[i] = &screenbuf[ (h - i - 1)*rowbytes ]; -#ifndef NOPNG writePng(outfn_conv.str().c_str(), rows, w, h); -#else // NOPNG - debMsgStd("ntlWorld::renderScene",DM_NOTIFY, "No PNG linked, no picture...", 1); -#endif // NOPNG } @@ -838,7 +848,7 @@ int ntlWorld::renderScene( void ) timeEnd = getTime(); char resout[1024]; - snprintf(resout,1024, "NTL Done %s, frame %d/%d (%s scene, %s raytracing, %s total, %d shades, %d i.s.'s)!\n", + snprintf(resout,1024, "NTL Done %s, frame %d/%d (took %s scene, %s raytracing, %s total, %d shades, %d i.s.'s)!\n", outfn_conv.str().c_str(), (glob->getAniCount()), (glob->getAniFrames()+1), getTimeString(totalStart-timeStart).c_str(), getTimeString(rendEnd-rendStart).c_str(), getTimeString(timeEnd-timeStart).c_str(), glob->getCounterShades(), diff --git a/intern/elbeem/intern/ntl_world.h b/intern/elbeem/intern/ntl_world.h index 18f388a5af4..9c5324cfe8f 100644 --- a/intern/elbeem/intern/ntl_world.h +++ b/intern/elbeem/intern/ntl_world.h @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Main renderer class * diff --git a/intern/elbeem/intern/parametrizer.cpp b/intern/elbeem/intern/parametrizer.cpp index b9dc1c36239..ae2bc7f3079 100644 --- a/intern/elbeem/intern/parametrizer.cpp +++ b/intern/elbeem/intern/parametrizer.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Parameter calculator for the LBM Solver class * @@ -57,7 +57,7 @@ Parametrizer::Parametrizer( void ) : mcAniFrameTime(0.0001), mTimeStepScale(1.0), mAniStart(0.0), - mExtent(1.0, 1.0, 1.0), //mSurfaceTension( 0.0 ), + //mExtent(1.0, 1.0, 1.0), //mSurfaceTension( 0.0 ), mDensity(1000.0), mGStar(0.0001), mFluidVolumeHeight(0.0), mSimulationMaxSpeed(0.0), mTadapMaxOmega(2.0), mTadapMaxSpeed(0.1), mTadapLevels(1), @@ -184,6 +184,7 @@ ParamFloat Parametrizer::calculateCellSize(void) int maxsize = mSizex; // get max size if(mSizey>maxsize) maxsize = mSizey; if(mSizez>maxsize) maxsize = mSizez; + maxsize = mSizez; // take along gravity dir for now! ParamFloat cellSize = 1.0 / (ParamFloat)maxsize; return cellSize; } @@ -252,9 +253,9 @@ int Parametrizer::calculateNoOfSteps( ParamFloat timelen ) { } /*! get extent of the domain = (1,1,1) if parametrizer not used, (x,y,z) [m] otherwise */ -ParamVec Parametrizer::calculateExtent( void ) { - return mExtent; -} +//ParamVec Parametrizer::calculateExtent( void ) { + //return mExtent; +//} /*! get (scaled) surface tension */ //ParamFloat Parametrizer::calculateSurfaceTension( void ) { return mSurfaceTension; } @@ -313,6 +314,7 @@ bool Parametrizer::calculateAllMissingValues( double time, bool silent ) int maxsize = mSizex; // get max size if(mSizey>maxsize) maxsize = mSizey; if(mSizez>maxsize) maxsize = mSizez; + maxsize = mSizez; // take along gravity dir for now! mCellSize = ( mDomainSize * calculateCellSize() ); // sets mCellSize if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," max domain resolution="<<(maxsize)<<" cells , cellsize="<0)) { // dump to binary file std::ostringstream boutfilename(""); - boutfilename << outfilename <<"_particles_" << frameNrStr<< ".gz"; + boutfilename << outfilename <<"_particles_" << frameNrStr; + if(glob_mpactive) { + if(glob_mpindex>0) { boutfilename << "mp"<getName() <<", particles:"<0.) { + debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"t="<0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) { // dump to binary file diff --git a/intern/elbeem/intern/particletracer.h b/intern/elbeem/intern/particletracer.h index f08bd58dc09..aae92aa1dea 100644 --- a/intern/elbeem/intern/particletracer.h +++ b/intern/elbeem/intern/particletracer.h @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Particle Viewer/Tracer * @@ -37,12 +37,14 @@ class ParticleObject public: //! Standard constructor inline ParticleObject(ntlVec3Gfx mp) : - mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0) { mId = ParticleObjectIdCnt++; }; + mPos(mp),mVel(0.0), mSize(1.0), mStatus(0),mLifeTime(0),mpNext(NULL) + { mId = ParticleObjectIdCnt++; }; //! Copy constructor inline ParticleObject(const ParticleObject &a) : mPos(a.mPos), mVel(a.mVel), mSize(a.mSize), mStatus(a.mStatus), - mLifeTime(a.mLifeTime) { mId = ParticleObjectIdCnt++; }; + mLifeTime(a.mLifeTime), mpNext(NULL) + { mId = ParticleObjectIdCnt++; }; //! Destructor inline ~ParticleObject() { /* empty */ }; @@ -70,6 +72,10 @@ class ParticleObject inline gfxReal getSize() { return mSize; } inline void setSize(gfxReal set) { mSize=set; } + //! get/set next pointer + inline ParticleObject* getNext() { return mpNext; } + inline void setNext(ParticleObject* set) { mpNext=set; } + //! get whole flags inline int getFlags() const { return mStatus; } //! get status (higher byte) @@ -101,6 +107,10 @@ class ParticleObject inline int getId() const { return mId; } + static inline float getMass(float size) { + return 4.0/3.0 * M_PI* (size)*(size)*(size); // mass: 4/3 pi r^3 rho + } + protected: /*! only for debugging */ @@ -115,6 +125,9 @@ class ParticleObject int mStatus; /*! count survived time steps */ float mLifeTime; + + /* for list constructions */ + ParticleObject *mpNext; }; @@ -130,6 +143,8 @@ class ParticleTracer : //! add a particle at this position void addParticle(float x, float y, float z); + //! add/copy a particle from inited struct + void addFullParticle(ParticleObject &np); //! draw the particle array void draw(); diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp index 3f954f1a4c7..fcec39810f9 100644 --- a/intern/elbeem/intern/simulation_object.cpp +++ b/intern/elbeem/intern/simulation_object.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Basic interface for all simulation modules * @@ -102,6 +102,7 @@ void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) { /****************************************************************************** * simluation interface: initialize simulation using the given configuration file *****************************************************************************/ +extern int glob_mpnum; int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) { if(! isSimworldOk() ) return 1; @@ -129,6 +130,8 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) } debMsgStd("SimulationObject::initialized",DM_MSG,"IdStr:"<getIdString() <<" LBM solver! ", 2); + mpParts = new ParticleTracer(); + // for non-param simulations mpLbm->setParametrizer( mpParam ); mpParam->setAttrList( getAttributeList() ); @@ -136,8 +139,8 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) mpParam->parseAttrList(); mpLbm->setAttrList( getAttributeList() ); + mpLbm->setSwsAttrList( getSwsAttributeList() ); mpLbm->parseAttrList(); - mpParts = new ParticleTracer(); mpParts->parseAttrList( getAttributeList() ); if(! isSimworldOk() ) return 3; @@ -163,6 +166,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) if(mpElbeemSettings->outputPath) this->mOutFilename = string(mpElbeemSettings->outputPath); mpLbm->initDomainTrafo( mpElbeemSettings->surfaceTrafo ); mpLbm->setSmoothing(1.0 * mpElbeemSettings->surfaceSmoothing, 1.0 * mpElbeemSettings->surfaceSmoothing); + mpLbm->setIsoSubdivs(mpElbeemSettings->surfaceSubdivs); mpLbm->setSizeX(mpElbeemSettings->resolutionxyz); mpLbm->setSizeY(mpElbeemSettings->resolutionxyz); mpLbm->setSizeZ(mpElbeemSettings->resolutionxyz); @@ -173,14 +177,14 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) mpParts->setNumInitialParticles(mpElbeemSettings->numTracerParticles); string dinitType = string("no"); - if (mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_PARTSLIP) dinitType = string("part"); - else if(mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_FREESLIP) dinitType = string("free"); - else /*if(mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_NOSLIP)*/ dinitType = string("no"); + if (mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_PARTSLIP) dinitType = string("part"); + else if(mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_FREESLIP) dinitType = string("free"); + else /*if(mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_NOSLIP)*/ dinitType = string("no"); mpLbm->setDomainBound(dinitType); - mpLbm->setDomainPartSlip(mpElbeemSettings->obstaclePartslip); + mpLbm->setDomainPartSlip(mpElbeemSettings->domainobsPartslip); mpLbm->setDumpVelocities(mpElbeemSettings->generateVertexVectors); mpLbm->setFarFieldSize(mpElbeemSettings->farFieldSize); - debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<generateVertexVectors<<","<getDumpVelocities(), 9 ); + debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<generateVertexVectors<<","<getDumpVelocities(), 9 ); debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<getGenerateParticles(),10); } @@ -193,63 +197,68 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob) if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverPostin status"); mPanic=true; return 15; } // print cell type stats - const int jmax = sizeof(CellFlagType)*8; - int totalCells = 0; - int flagCount[jmax]; - for(int j=0; jgetFirstCell(); - for(; mpLbm->noEndCell( cid ); - mpLbm->advanceCell( cid ) ) { - int flag = mpLbm->getCellFlag(cid,0); - int flag2 = mpLbm->getCellFlag(cid,1); - if(flag != flag2) { - diffInits++; + bool printStats = true; + if(glob_mpnum>0) printStats=false; // skip in this case + if(printStats) { + const int jmax = sizeof(CellFlagType)*8; + int totalCells = 0; + int flagCount[jmax]; + for(int j=0; jgetFirstCell(); + for(; mpLbm->noEndCell( cid ); + mpLbm->advanceCell( cid ) ) { + int flag = mpLbm->getCellFlag(cid,0); + int flag2 = mpLbm->getCellFlag(cid,1); + if(flag != flag2) { + diffInits++; + } + for(int j=0; jdeleteCellIterator( &cid ); + mpLbm->deleteCellIterator( &cid ); - char charNl = '\n'; - debugOutNnl("SimulationObject::initializeLbmSimulation celltype stats: " <0) { - out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1<0) { + out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1< 0) { + debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="< 0) { - debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="<setStart( mGeoStart ); - mpParts->setEnd( mGeoEnd ); + //mpParts->setStart( mGeoStart );? mpParts->setEnd( mGeoEnd );? + mpParts->setStart( mpLbm->getGeoStart() ); + mpParts->setEnd( mpLbm->getGeoEnd() ); mpParts->setCastShadows( false ); mpParts->setReceiveShadows( false ); mpParts->searchMaterial( glob->getMaterials() ); // this has to be inited here - before, the values might be unknown - ntlGeometryObject *surf = mpLbm->getSurfaceGeoObj(); + IsoSurface *surf = mpLbm->getSurfaceGeoObj(); if(surf) { surf->setName( "final" ); // final surface mesh // warning - this might cause overwriting effects for multiple sims and geom dump... diff --git a/intern/elbeem/intern/simulation_object.h b/intern/elbeem/intern/simulation_object.h index 200cad85221..56d7f20e7cd 100644 --- a/intern/elbeem/intern/simulation_object.h +++ b/intern/elbeem/intern/simulation_object.h @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Basic interface for all simulation modules * diff --git a/intern/elbeem/intern/solver_adap.cpp b/intern/elbeem/intern/solver_adap.cpp index 99892081fca..ef516a578bd 100644 --- a/intern/elbeem/intern/solver_adap.cpp +++ b/intern/elbeem/intern/solver_adap.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004,2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Adaptivity functions * @@ -925,6 +925,11 @@ void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int d +// required globals +extern bool glob_mpactive; +extern int glob_mpnum, glob_mpindex; +#define MPTADAP_INTERV 4 + /*****************************************************************************/ /*! change the size of the LBM time step */ /*****************************************************************************/ @@ -934,7 +939,7 @@ void LbmFsgrSolver::adaptTimestep() { bool rescale = false; // do any rescale at all? LbmFloat scaleFac = -1.0; // timestep scaling - if(this->mPanic) return; + if(mPanic) return; LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS]; LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS]; @@ -942,62 +947,71 @@ void LbmFsgrSolver::adaptTimestep() { levOldOmega[lev] = mLevel[lev].omega; levOldStepsize[lev] = mLevel[lev].timestep; } - //if(mTimeSwitchCounts>0){ errMsg("DEB CSKIP",""); return; } // DEBUG const LbmFloat reduceFac = 0.8; // modify time step by 20%, TODO? do multiple times for large changes? LbmFloat diffPercent = 0.05; // dont scale if less than 5% - LbmFloat allowMax = this->mpParam->getTadapMaxSpeed(); // maximum allowed velocity - LbmFloat nextmax = this->mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity); + LbmFloat allowMax = mpParam->getTadapMaxSpeed(); // maximum allowed velocity + LbmFloat nextmax = mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity); - //newdt = this->mpParam->getTimestep() * (allowMax/nextmax); - LbmFloat newdt = this->mpParam->getTimestep(); // newtr + // sync nextmax +#if LBM_INCLUDE_TESTSOLVERS==1 + if(glob_mpactive) { + if(mLevel[mMaxRefine].lsteps % MPTADAP_INTERV != MPTADAP_INTERV-1) { + debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "mpact:"<getTimestep(); // newtr if(nextmax > allowMax/reduceFac) { mTimeMaxvelStepCnt++; } else { mTimeMaxvelStepCnt=0; } // emergency, or 10 steps with high vel if((mTimeMaxvelStepCnt>5) || (nextmax> (1.0/3.0)) || (mForceTimeStepReduce) ) { - //if(nextmax > allowMax/reduceFac) { - newdt = this->mpParam->getTimestep() * reduceFac; + newdt = mpParam->getTimestep() * reduceFac; } else { if(nextmaxmpParam->getTimestep() / reduceFac; + newdt = mpParam->getTimestep() / reduceFac; } } // newtr - //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<mpParam->getMaxTimestep()<<" min"<mpParam->getMinTimestep()<<" diff"<getMaxTimestep()<<" min"<getMinTimestep()<<" diff"<mpParam->getTimestep() ) // DEBUG + //if( newdt < mpParam->getTimestep() ) // DEBUG LbmFloat rhoAvg = mCurrentMass/mCurrentVolume; - if( (newdt<=this->mpParam->getMaxTimestep()) && (newdt>=this->mpParam->getMinTimestep()) - && (dtdiff>(this->mpParam->getTimestep()*diffPercent)) ) { + if( (newdt<=mpParam->getMaxTimestep()) && (newdt>=mpParam->getMinTimestep()) + && (dtdiff>(mpParam->getTimestep()*diffPercent)) ) { if((newdt>levOldStepsize[mMaxRefine])&&(mTimestepReduceLock)) { // wait some more... //debMsgNnl("LbmFsgrSolver::TAdp",DM_NOTIFY," Delayed... "<mpParam->setDesiredTimestep( newdt ); + mpParam->setDesiredTimestep( newdt ); rescale = true; - if(!this->mSilent) { + if(!mSilent) { debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10); - debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<mpParam->getSimulationMaxSpeed()<<" next:"<mStepCnt, 10 ); - debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<< - "rhoAvg="<getSimulationMaxSpeed()<<" next:"<mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) || - ((this->mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){ + if( ((mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) || + ((mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){ rescale = true; minCutoff = false; - newdt = tadtogSwitch * this->mpParam->getTimestep(); - this->mpParam->setDesiredTimestep( newdt ); + newdt = tadtogSwitch * mpParam->getTimestep(); + mpParam->setDesiredTimestep( newdt ); } else - if( ((this->mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) || - ((this->mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){ + if( ((mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) || + ((mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){ rescale = true; minCutoff = false; - newdt = this->mpParam->getTimestep()/tadtogSwitch ; - this->mpParam->setDesiredTimestep( newdt ); + newdt = mpParam->getTimestep()/tadtogSwitch ; + mpParam->setDesiredTimestep( newdt ); } else { rescale = false; minCutoff = false; } @@ -1027,23 +1041,25 @@ void LbmFsgrSolver::adaptTimestep() { // test mass rescale - scaleFac = newdt/this->mpParam->getTimestep(); + scaleFac = newdt/mpParam->getTimestep(); if(rescale) { // perform acutal rescaling... mTimeMaxvelStepCnt=0; mForceTimeStepReduce = false; // FIXME - approximate by averaging, take gravity direction here? - mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3; + //mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3; + // use z as gravity direction + mTimestepReduceLock = 4*mLevel[mMaxRefine].lSizez; mTimeSwitchCounts++; - this->mpParam->calculateAllMissingValues( mSimulationTime, this->mSilent ); + mpParam->calculateAllMissingValues( mSimulationTime, mSilent ); recalculateObjectSpeeds(); // calc omega, force for all levels mLastOmega=1e10; mLastGravity=1e10; initLevelOmegas(); - if(this->mpParam->getTimestep()mpParam->getTimestep(); - if(this->mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = this->mpParam->getTimestep(); + if(mpParam->getTimestep()getTimestep(); + if(mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = mpParam->getTimestep(); // this might be called during init, before we have any particles if(mpParticles) { mpParticles->adaptPartTimestep(scaleFac); } @@ -1058,8 +1074,8 @@ void LbmFsgrSolver::adaptTimestep() { LbmFloat newSteptime = mLevel[lev].timestep; LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]); - if(!this->mSilent) { - debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<mStepCnt<< + if(!mSilent) { + debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<mOmega<<" gStar:"<mpParam->getCurrentGStar() , 10); + debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep changed by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<getCurrentGStar()<<", step:"<mName<<" step:"<mStepCnt<<" newdt="<=0 ; lev--) { @@ -1220,12 +1236,12 @@ void LbmFsgrSolver::adaptTimestep() { // collide on current set LbmFloat rho, ux,uy,uz; rho=0.0; ux = uy = uz = 0.0; - for(int l=0; lcDfNum; l++) { + for(int l=0; ldfDvecX[l]*m); - uy += (this->dfDvecY[l]*m); - uz += (this->dfDvecZ[l]*m); + ux += (dfDvecX[l]*m); + uy += (dfDvecY[l]*m); + uz += (dfDvecZ[l]*m); } #ifndef WIN32 if (!finite(rho)) { @@ -1242,8 +1258,8 @@ void LbmFsgrSolver::adaptTimestep() { ux *= cfac; uy *= cfac; uz *= cfac; - for(int l=0; lcDfNum; l++) { - QCELL(lev, i, j, k, workSet, l) = this->getCollideEq(l, rho, ux,uy,uz); } + for(int l=0; lLbmCell LbmCell; - protected: //! internal quick print function (for debugging) @@ -316,6 +322,11 @@ class LbmFsgrSolver : void reinitFlags( int workSet ); //! mass dist weights LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l); + //! compute surface normals: fluid, fluid more accurate, and for obstacles + void computeFluidSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); + void computeFluidSurfaceNormalAcc(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); + void computeObstacleSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); + void computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret); //! add point to mListNewInter list LBM_INLINED void addToNewInterList( int ni, int nj, int nk ); //! cell is interpolated from coarse level (inited into set, source sets are determined by t) @@ -327,6 +338,8 @@ class LbmFsgrSolver : LBM_INLINED int getForZMin1(); LBM_INLINED int getForZMaxBnd(int lev); LBM_INLINED int getForZMax1(int lev); + LBM_INLINED bool checkDomainBounds(int lev,int i,int j,int k); + LBM_INLINED bool checkDomainBoundsPos(int lev,LbmVec pos); // touch grid and flags once void preinitGrids(); @@ -361,8 +374,6 @@ class LbmFsgrSolver : LbmFloat mFVArea; bool mUpdateFVHeight; - //! require some geo setup from the viz? - //int mGfxGeoSetup; //! force quit for gfx LbmFloat mGfxEndTime; //! smoother surface initialization? @@ -500,10 +511,21 @@ class LbmFsgrSolver : void handleCpdata(); void cpDebugDisplay(int dispset); - void testXchng(); + int mMpNum,mMpIndex; + int mOrgSizeX; + LbmFloat mOrgStartX; + LbmFloat mOrgEndX; + void mrSetup(); + void mrExchange(); + void mrIsoExchange(); + LbmFloat mrInitTadap(LbmFloat max); + void gcFillBuffer( LbmGridConnector *gc, int *retSizeCnt, const int *bdfs); + void gcUnpackBuffer(LbmGridConnector *gc, int *retSizeCnt, const int *bdfs); public: // needed for testdata void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz); + // mptest + int getMpIndex() { return mMpIndex; }; # endif // LBM_INCLUDE_TESTSOLVERS==1 // former LbmModelLBGK functions @@ -615,7 +637,7 @@ class LbmFsgrSolver : STCON LbmFloat dfLength[ 19 ]; /*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */ - static LbmFloat dfEquil[ 19 ]; + static LbmFloat dfEquil[ dTotalNum ]; /*! arrays for les model coefficients */ static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ]; @@ -701,7 +723,7 @@ class LbmFsgrSolver : STCON LbmFloat dfLength[ 9 ]; /* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */ - static LbmFloat dfEquil[ 9 ]; + static LbmFloat dfEquil[ dTotalNum ]; /*! arrays for les model coefficients */ static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ]; @@ -784,7 +806,8 @@ class LbmFsgrSolver : // general defines ----------------------------------------------------------------------------------------------- -#define TESTFLAG(flag, compflag) ((flag & compflag)==compflag) +// replace TESTFLAG +#define FLAGISEXACT(flag, compflag) ((flag & compflag)==compflag) #if LBMDIM==2 #define dC 0 @@ -943,6 +966,41 @@ int LbmFsgrSolver::getForZMax1(int lev) { return mLevel[lev].lSizez -1; } +bool LbmFsgrSolver::checkDomainBounds(int lev,int i,int j,int k) { + if(i<0) return false; + if(j<0) return false; + if(k<0) return false; + if(i>mLevel[lev].lSizex-1) return false; + if(j>mLevel[lev].lSizey-1) return false; + if(k>mLevel[lev].lSizez-1) return false; + return true; +} +bool LbmFsgrSolver::checkDomainBoundsPos(int lev,LbmVec pos) { + const int i= (int)pos[0]; + if(i<0) return false; + if(i>mLevel[lev].lSizex-1) return false; + const int j= (int)pos[1]; + if(j<0) return false; + if(j>mLevel[lev].lSizey-1) return false; + const int k= (int)pos[2]; + if(k<0) return false; + if(k>mLevel[lev].lSizez-1) return false; + return true; +} + +void LbmFsgrSolver::initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass) { + LbmFloat *ccel = &QCELL(level ,i,j,k, workSet,0); + LbmFloat nrho = 0.0; + FORDF0 { nrho += RAC(ccel,l); } + if(initMass) { + RAC(ccel,dMass) = nrho; + RAC(ccel, dFfrac) = 1.; + } else { + // preinited, e.g. from reinitFlags + RAC(ccel, dFfrac) = RAC(ccel, dMass)/nrho; + RAC(ccel, dFlux) = FLUX_INIT; + } +} #endif diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp index 1a19e1deabf..6f891310224 100644 --- a/intern/elbeem/intern/solver_init.cpp +++ b/intern/elbeem/intern/solver_init.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004,2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Standard LBM Factory implementation * @@ -13,9 +13,6 @@ // for geo init FGI_ defines #include "elbeem.h" -// all boundary types at once -#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW ) - // helper for 2d init #define SWAPYZ(vec) { \ const LbmFloat tmp = (vec)[2]; \ @@ -173,7 +170,7 @@ }; /* precalculated equilibrium dfs, inited in lbmsolver constructor */ - LbmFloat LbmFsgrSolver::dfEquil[ cDfNum ]; + LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ]; #else // end LBMDIM==3 , LBMDIM==2 @@ -290,12 +287,17 @@ }; /* precalculated equilibrium dfs, inited in lbmsolver constructor */ - LbmFloat LbmFsgrSolver::dfEquil[ cDfNum ]; + LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ]; // D2Q9 end #endif // LBMDIM==2 +// required globals +extern bool glob_mpactive; +extern int glob_mpnum, glob_mpindex; + + /****************************************************************************** * Lbm Constructor *****************************************************************************/ @@ -318,7 +320,7 @@ LbmFsgrSolver::LbmFsgrSolver() : mIsoWeightMethod(1), mMaxRefine(1), mDfScaleUp(-1.0), mDfScaleDown(-1.0), - mInitialCsmago(0.03), // set to 0.02 for mMaxRefine==0 below + mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03 mDebugOmegaRet(0.0), mLastOmega(1e10), mLastGravity(1e10), mNumInvIfTotal(0), mNumFsgrChanges(0), @@ -329,38 +331,45 @@ LbmFsgrSolver::LbmFsgrSolver() : // not much to do here... #if LBM_INCLUDE_TESTSOLVERS==1 mpTest = new LbmTestdata(); -#endif // ELBEEM_PLUGIN!=1 - this->mpIso = new IsoSurface( this->mIsoValue ); + mMpNum = mMpIndex = 0; + mOrgSizeX = 0; + mOrgStartX = 0.; + mOrgEndX = 0.; +#endif // LBM_INCLUDE_TESTSOLVERS!=1 + mpIso = new IsoSurface( mIsoValue ); // init equilibrium dist. func LbmFloat rho=1.0; FORDF0 { - this->dfEquil[l] = this->getCollideEq( l,rho, 0.0, 0.0, 0.0); + dfEquil[l] = this->getCollideEq( l,rho, 0.0, 0.0, 0.0); } + dfEquil[dMass] = 1.; + dfEquil[dFfrac] = 1.; + dfEquil[dFlux] = FLUX_INIT; // init LES int odm = 0; for(int m=0; mcDfNum; l++) { + for(int l=0; llesCoeffDiag[m][l] = this->lesCoeffOffdiag[m][l] = 0.0; } } for(int m=0; mcDfNum; l++) { + for(int l=1; ldfDvecX[l]; break; - case 1: em = this->dfDvecY[l]; break; - case 2: em = this->dfDvecZ[l]; break; + case 0: em = dfDvecX[l]; break; + case 1: em = dfDvecY[l]; break; + case 2: em = dfDvecZ[l]; break; default: em = -1.0; errFatal("SMAGO1","err m="<dfDvecX[l]; break; - case 1: en = this->dfDvecY[l]; break; - case 2: en = this->dfDvecZ[l]; break; + case 0: en = dfDvecX[l]; break; + case 1: en = dfDvecY[l]; break; + case 2: en = dfDvecZ[l]; break; default: en = -1.0; errFatal("SMAGO2","err n="<dfDvecX[this->dfInv[l]], this->dfDvecY[this->dfInv[l]], this->dfDvecZ[this->dfInv[l]] ) + LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] ) ) * -1.0; } @@ -395,15 +404,15 @@ LbmFsgrSolver::LbmFsgrSolver() : #if ELBEEM_PLUGIN!=1 errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!"); #endif - for(int n=0;(ncDirNum); n++) { mGaussw[n] = 0.0; } - //for(int n=0;(ncDirNum); n++) { - for(int n=0;(ncDfNum); n++) { - const LbmFloat d = norm(LbmVec(this->dfVecX[n], this->dfVecY[n], this->dfVecZ[n])); + for(int n=0;(ncDirNum); n++) { + for(int n=0;(nmInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; } + if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; } #if COMPRESSGRIDS==1 delete [] mLevel[mMaxRefine].mprsCells[1]; mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL; @@ -426,16 +435,13 @@ LbmFsgrSolver::~LbmFsgrSolver() if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s]; } } - delete this->mpIso; + delete mpIso; if(mpPreviewSurface) delete mpPreviewSurface; - -#if LBM_INCLUDE_TESTSOLVERS==1 // cleanup done during scene deletion... -#endif // ELBEEM_PLUGIN!=1 // always output performance estimate debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5); - if(!this->mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10); + if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10); } @@ -449,19 +455,19 @@ void LbmFsgrSolver::parseAttrList() LbmSolverInterface::parseStdAttrList(); string matIso("default"); - matIso = this->mpAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false ); - this->mpIso->setMaterialName( matIso ); - this->mOutputSurfacePreview = this->mpAttrs->readInt("surfacepreview", this->mOutputSurfacePreview, "SimulationLbm","this->mOutputSurfacePreview", false ); - mTimeAdap = this->mpAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false ); - this->mDomainBound = this->mpAttrs->readString("domainbound", this->mDomainBound, "SimulationLbm","mDomainBound", false ); - this->mDomainPartSlipValue = this->mpAttrs->readFloat("domainpartslip", this->mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false ); + matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false ); + mpIso->setMaterialName( matIso ); + mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false ); + mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false ); + mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false ); + mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false ); - mIsoWeightMethod= this->mpAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false ); - mInitSurfaceSmoothing = this->mpAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false ); - this->mSmoothSurface = this->mpAttrs->readFloat("smoothsurface", this->mSmoothSurface, "SimulationLbm","mSmoothSurface", false ); - this->mSmoothNormals = this->mpAttrs->readFloat("smoothnormals", this->mSmoothNormals, "SimulationLbm","mSmoothNormals", false ); + mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false ); + mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false ); + mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false ); + mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false ); - mFsSurfGenSetting = this->mpAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false ); + mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false ); if(mFsSurfGenSetting==-1) { // all on mFsSurfGenSetting = @@ -470,46 +476,61 @@ void LbmFsgrSolver::parseAttrList() } // refinement - mMaxRefine = this->mRefinementDesired; - mMaxRefine = this->mpAttrs->readInt("maxrefine", mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false); + mMaxRefine = mRefinementDesired; + mMaxRefine = mpSifAttrs->readInt("maxrefine", mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false); if(mMaxRefine<0) mMaxRefine=0; if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1; - mDisableStandingFluidInit = this->mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false); - mInit2dYZ = this->mpAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false); - mForceTadapRefine = this->mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false); + mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false); + mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false); + mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false); // demo mode settings - mFVHeight = this->mpAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false ); + mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false ); // FIXME check needed? - mFVArea = this->mpAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false ); + mFVArea = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false ); // debugging - skip some time... double starttimeskip = 0.; - starttimeskip = this->mpAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false ); + starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false ); mSimulationTime += starttimeskip; if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<parseTestdataAttrList(this->mpAttrs); + mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false); + mpTest->parseTestdataAttrList(mpSifAttrs); #ifdef ELBEEM_PLUGIN mUseTestdata=1; // DEBUG #endif // ELBEEM_PLUGIN - errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<readInt("mpnum", mMpNum ,"LbmFsgrSolver", "mMpNum", false); + mMpIndex = mpSifAttrs->readInt("mpindex", mMpIndex ,"LbmFsgrSolver", "mMpIndex", false); + if(glob_mpactive) { + // used instead... + mMpNum = glob_mpnum; + mMpIndex = glob_mpindex; + } else { + glob_mpnum = mMpNum; + glob_mpindex = 0; + } + errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<0) { + mUseTestdata=1; // needed in this case... + } + + errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<=2.) mUseTestdata=1; // equiv. to test solver check #endif // LBM_INCLUDE_TESTSOLVERS!=1 //if(mUseTestdata) { mMaxRefine=0; } // force fsgr off if(mMaxRefine==0) mInitialCsmago=0.02; - mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false ); + mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false ); // deprecated! float mInitialCsmagoCoarse = 0.0; - mInitialCsmagoCoarse = this->mpAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false ); + mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false ); #if USE_LES==1 #else // USE_LES==1 debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2); @@ -524,15 +545,15 @@ void LbmFsgrSolver::parseAttrList() void LbmFsgrSolver::initLevelOmegas() { // no explicit settings - this->mOmega = this->mpParam->calculateOmega(mSimulationTime); - this->mGravity = vec2L( this->mpParam->calculateGravity(mSimulationTime) ); - this->mSurfaceTension = 0.; //this->mpParam->calculateSurfaceTension(); // unused + mOmega = mpParam->calculateOmega(mSimulationTime); + mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) ); + mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused if(mInit2dYZ) { SWAPYZ(mGravity); } // check if last init was ok - LbmFloat gravDelta = norm(this->mGravity-mLastGravity); - //errMsg("ChannelAnimDebug","t:"<mOmega<<" - lom:"<mGravity<<" - "<mOmega == mLastOmega) && (gravDelta<=0.0)) return; + LbmFloat gravDelta = norm(mGravity-mLastGravity); + //errMsg("ChannelAnimDebug","t:"<mOmega; - mLevel[i].timestep = this->mpParam->getTimestep(); + mLevel[i].omega = mOmega; + mLevel[i].timestep = mpParam->getTimestep(); mLevel[i].lcsmago = mInitialCsmago; //CSMAGO_INITIAL; mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago; mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0); } // init all sub levels + LbmFloat coarseCsmago = mInitialCsmago; + coarseCsmago = 0.04; // try stabilizing for(int i=mMaxRefine-1; i>=0; i--) { //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5; double nomega = 0.5 * ( (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5; nomega = 1.0/nomega; mLevel[i].omega = (LbmFloat)nomega; mLevel[i].timestep = 2.0 * mLevel[i+1].timestep; - mLevel[i].lcsmago = mInitialCsmago; + mLevel[i].lcsmago = coarseCsmago; mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago; mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0); } // for lbgk - mLevel[ mMaxRefine ].gravity = this->mGravity / mLevel[ mMaxRefine ].omega; + mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega; for(int i=mMaxRefine-1; i>=0; i--) { // should be the same on all levels... // for lbgk mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega; } - mLastOmega = this->mOmega; - mLastGravity = this->mGravity; + mLastOmega = mOmega; + mLastGravity = mGravity; // debug? invalidate old values... - this->mGravity = -100.0; - this->mOmega = -100.0; + mGravity = -100.0; + mOmega = -100.0; for(int i=0; i<=mMaxRefine; i++) { - if(!this->mSilent) { + if(!mSilent) { errMsg("LbmFsgrSolver", "Level init "<mInitDone) { + if(!mInitDone) { debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<mInitDone<<" "<<(void*)this,1); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<0) { - this->mSizex *= mCppfStage; - this->mSizey *= mCppfStage; - this->mSizez *= mCppfStage; + mSizex *= mCppfStage; + mSizey *= mCppfStage; + mSizez *= mCppfStage; } // size inits to force cubic cells and mult4 level dimensions // and make sure we dont allocate too much... bool memOk = false; - int orgSx = this->mSizex; - int orgSy = this->mSizey; - int orgSz = this->mSizez; + int orgSx = mSizex; + int orgSy = mSizey; + int orgSz = mSizez; double sizeReduction = 1.0; double memEstFromFunc = -1.0; string memreqStr(""); + bool firstMInit = true; + int minitTries=0; while(!memOk) { - initGridSizes( this->mSizex, this->mSizey, this->mSizez, - this->mvGeoStart, this->mvGeoEnd, mMaxRefine, PARALLEL); - calculateMemreqEstimate( this->mSizex, this->mSizey, this->mSizez, mMaxRefine, &memEstFromFunc, &memreqStr ); + minitTries++; + initGridSizes( mSizex, mSizey, mSizez, + mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL); + + // MPT +#if LBM_INCLUDE_TESTSOLVERS==1 + if(firstMInit) { + mrSetup(); + } +#endif // LBM_INCLUDE_TESTSOLVERS==1 + firstMInit=false; + + calculateMemreqEstimate( mSizex, mSizey, mSizez, + mMaxRefine, mFarFieldSize, &memEstFromFunc, &memreqStr ); double memLimit; + string memLimStr("-"); if(sizeof(void*)==4) { // 32bit system, limit to 2GB memLimit = 2.0* 1024.0*1024.0*1024.0; + memLimStr = string("2GB"); } else { // 64bit, just take 16GB as limit for now... memLimit = 16.0* 1024.0*1024.0*1024.0; + memLimStr = string("16GB"); } if(memEstFromFunc>memLimit) { sizeReduction *= 0.9; - this->mSizex = (int)(orgSx * sizeReduction); - this->mSizey = (int)(orgSy * sizeReduction); - this->mSizez = (int)(orgSz * sizeReduction); - debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<mSizex<<" Y:"<mSizey<<" Z:"<mSizez, 3 ); + mSizex = (int)(orgSx * sizeReduction); + mSizey = (int)(orgSy * sizeReduction); + mSizez = (int)(orgSz * sizeReduction); + debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<< + //memEstFromFunc<<"/"<mPreviewFactor = (LbmFloat)this->mOutputSurfacePreview / (LbmFloat)this->mSizex; - debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<mSizex<<" Y:"<mSizey<<" Z:"<mSizez<< - ", Domain: "<mvGeoStart<<":"<mvGeoEnd<<", "<<(this->mvGeoEnd-this->mvGeoStart)<< - ", Intbytes: "<< sizeof(void*) << + mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex; + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<mpParam->setSize(this->mSizex, this->mSizey, this->mSizez); + mpParam->setSize(mSizex, mSizey, mSizez); + if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); } debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: " <<"LBM_EPSILON="<mSizez = 1; + if(LBMDIM == 2) mSizez = 1; - this->mpParam->setSimulationMaxSpeed(0.0); - if(mFVHeight>0.0) this->mpParam->setFluidVolumeHeight(mFVHeight); - this->mpParam->setTadapLevels( mMaxRefine+1 ); + mpParam->setSimulationMaxSpeed(0.0); + if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight); + mpParam->setTadapLevels( mMaxRefine+1 ); if(mForceTadapRefine>mMaxRefine) { - this->mpParam->setTadapLevels( mForceTadapRefine+1 ); + mpParam->setTadapLevels( mForceTadapRefine+1 ); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<mpParam->calculateAllMissingValues(mSimulationTime, false)) { + if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) { errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR); return false; } @@ -706,9 +747,9 @@ bool LbmFsgrSolver::initializeSolverMemory() } // init sizes - mLevel[mMaxRefine].lSizex = this->mSizex; - mLevel[mMaxRefine].lSizey = this->mSizey; - mLevel[mMaxRefine].lSizez = this->mSizez; + mLevel[mMaxRefine].lSizex = mSizex; + mLevel[mMaxRefine].lSizey = mSizey; + mLevel[mMaxRefine].lSizez = mSizez; for(int i=mMaxRefine-1; i>=0; i--) { mLevel[i].lSizex = mLevel[i+1].lSizex/2; mLevel[i].lSizey = mLevel[i+1].lSizey/2; @@ -722,8 +763,8 @@ bool LbmFsgrSolver::initializeSolverMemory() } double ownMemCheck = 0.0; - mLevel[ mMaxRefine ].nodeSize = ((this->mvGeoEnd[0]-this->mvGeoStart[0]) / (LbmFloat)(this->mSizex)); - mLevel[ mMaxRefine ].simCellSize = this->mpParam->getCellSize(); + mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex)); + mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize(); mLevel[ mMaxRefine ].lcellfactor = 1.0; LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum); // +4 for safety ? @@ -759,8 +800,14 @@ bool LbmFsgrSolver::initializeSolverMemory() ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4); } - // isosurface memory - ownMemCheck += (3*sizeof(int)+sizeof(float)) * ((this->mSizex+2)*(this->mSizey+2)*(this->mSizez+2)); + // isosurface memory, use orig res values + if(mFarFieldSize>0.) { + ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) ); + } else { + // ignore 3 int slices... + ownMemCheck += (double)( ( sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) ); + } + // sanity check #if ELBEEM_PLUGIN!=1 if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) { @@ -782,31 +829,31 @@ bool LbmFsgrSolver::initializeSolverMemory() // calc omega, force for all levels initLevelOmegas(); - mMinTimestep = this->mpParam->getTimestep(); - mMaxTimestep = this->mpParam->getTimestep(); + mMinTimestep = mpParam->getTimestep(); + mMaxTimestep = mpParam->getTimestep(); // init isosurf - this->mpIso->setIsolevel( this->mIsoValue ); + mpIso->setIsolevel( 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->setMaterialName( mpIso->getMaterialName() ); + delete mpIso; + mpIso = mpTest; + if(mpTest->mFarfMode>0) { // 3d off mpTest->setIsolevel(-100.0); } else { - mpTest->setIsolevel( this->mIsoValue ); + mpTest->setIsolevel( mIsoValue ); } } -#endif // ELBEEM_PLUGIN!=1 +#endif // LBM_INCLUDE_TESTSOLVERS!=1 // 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 ); + if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.; + if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.; + mpIso->setSmoothSurface( mSmoothSurface * featureSize ); + mpIso->setSmoothNormals( mSmoothNormals * featureSize ); // init iso weight values mIsoWeightMethod int wcnt = 0; @@ -835,34 +882,54 @@ bool LbmFsgrSolver::initializeSolverMemory() } for(int i=0; i<27; i++) mIsoWeight[i] /= totw; - LbmVec isostart = vec2L(this->mvGeoStart); - LbmVec isoend = vec2L(this->mvGeoEnd); + LbmVec isostart = vec2L(mvGeoStart); + LbmVec isoend = vec2L(mvGeoEnd); int twodOff = 0; // 2d slices if(LBMDIM==2) { LbmFloat sn,se; - sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5; - se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5; + sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5; + se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5; isostart[2] = sn; isoend[2] = se; twodOff = 2; } - int isosx = this->mSizex+2; - int isosy = this->mSizey+2; - int isosz = this->mSizez+2+twodOff; + int isosx = mSizex+2; + int isosy = mSizey+2; + int isosz = mSizez+2+twodOff; // MPT #if LBM_INCLUDE_TESTSOLVERS==1 - if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) { - int xfac = 2; - isosx *= xfac; - isoend[0] = isostart[0] + (isoend[0]-isostart[0])*(LbmFloat)xfac; + //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) { + if( (mMpNum>0) && (mMpIndex==0) ) { + //? mpindex==0 + // restore original value for node0 + isosx = mOrgSizeX + 2; + isostart[0] = mOrgStartX; + isoend[0] = mOrgEndX; } + errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst" + << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize" + << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" "); + errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst" + << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize" + << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" "); #endif // LBM_INCLUDE_TESTSOLVERS==1 - //errMsg(" SETISO ", " "<mSizex+1.0))*0.5)<<" "<<(LbmFloat)(this->mSizex+1.0)<<" " ); + errMsg(" SETISO ", "iso "<setStart( vec2G(isostart) ); mpIso->setEnd( vec2G(isoend) ); LbmVec isodist = isoend-isostart; + + int isosubs = mIsoSubdivs; + if(mFarFieldSize>1.) { + errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!"); + isosubs = 1; + mpIso->setUseFulledgeArrays(true); + } + mpIso->setSubdivs(isosubs); + mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) ); // reset iso field @@ -876,7 +943,7 @@ bool LbmFsgrSolver::initializeSolverMemory() for(int lev=0; lev<=mMaxRefine; lev++) { FSGR_FORIJK_BOUNDS(lev) { RFLAG(lev,i,j,k,0) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage - if(!this->mAllfluid) { + if(!mAllfluid) { initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); } else { initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); @@ -884,38 +951,43 @@ bool LbmFsgrSolver::initializeSolverMemory() } } + if(LBMDIM==2) { - if(this->mOutputSurfacePreview) { + if(mOutputSurfacePreview) { errMsg("LbmFsgrSolver::init","No preview in 2D allowed!"); - this->mOutputSurfacePreview = 0; } + mOutputSurfacePreview = 0; } } + if((glob_mpactive) && (glob_mpindex>0)) { + mOutputSurfacePreview = 0; + } + #if LBM_USE_GUI==1 - if(this->mOutputSurfacePreview) { + if(mOutputSurfacePreview) { errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0"); - this->mOutputSurfacePreview = 0; } + mOutputSurfacePreview = 0; } #endif // LBM_USE_GUI==1 - if(this->mOutputSurfacePreview) { + if(mOutputSurfacePreview) { // same as normal one, but use reduced size - mpPreviewSurface = new IsoSurface( this->mIsoValue ); + mpPreviewSurface = new IsoSurface( mIsoValue ); mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() ); - mpPreviewSurface->setIsolevel( this->mIsoValue ); + mpPreviewSurface->setIsolevel( mIsoValue ); // usually dont display for rendering mpPreviewSurface->setVisible( false ); mpPreviewSurface->setStart( vec2G(isostart) ); mpPreviewSurface->setEnd( vec2G(isoend) ); LbmVec pisodist = isoend-isostart; - LbmFloat pfac = this->mPreviewFactor; - mpPreviewSurface->initializeIsosurface( (int)(pfac*this->mSizex)+2, (int)(pfac*this->mSizey)+2, (int)(pfac*this->mSizez)+2, vec2G(pisodist) ); - //mpPreviewSurface->setName( this->getName() + "preview" ); + LbmFloat pfac = mPreviewFactor; + mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) ); + //mpPreviewSurface->setName( getName() + "preview" ); mpPreviewSurface->setName( "preview" ); - debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*this->mSizex)<<","<<(pfac*this->mSizey)<<","<<(pfac*this->mSizez)<<" enabled",10); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10); } // init defaults mAvgNumUsedCells = 0; - this->mFixMass= 0.0; + mFixMass= 0.0; return true; } @@ -933,23 +1005,23 @@ bool LbmFsgrSolver::initializeSolverGrids() { CellFlagType domainBoundType = CFInvalid; // TODO use normal object types instad... - if(this->mDomainBound.find(string("free")) != string::npos) { + if(mDomainBound.find(string("free")) != string::npos) { domainBoundType = CFBnd | CFBndFreeslip; - debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<mDomainBound,10); - } else if(this->mDomainBound.find(string("part")) != string::npos) { + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<mDomainPartSlipValue<<"), value:"<mDomainBound,10); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<mDomainBound,10); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<mpGiObjects->size()); + int domainobj = (int)(mpGiObjects->size()); domainBoundType |= (domainobj<<24); //for(int i=0; i<(int)(domainobj+0); i++) { - //errMsg("GEOIN","i"<mpGiObjects)[i]->getName()); - //if((*this->mpGiObjects)[i] == this->mpIso) { //check... + //errMsg("GEOIN","i"<getName()); + //if((*mpGiObjects)[i] == mpIso) { //check... //} //} //errMsg("GEOIN"," dm "<<(domainBoundType>>24)); @@ -1050,19 +1122,19 @@ bool LbmFsgrSolver::initializeSolverGrids() { mInitialMass = 0.0; int inmCellCnt = 0; FSGR_FORIJK1(mMaxRefine) { - if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid) ) { + if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) { LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); } mInitialMass += fluidRho; inmCellCnt ++; - } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) { + } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) { mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass); inmCellCnt ++; } } mCurrentVolume = mCurrentMass = mInitialMass; - ParamVec cspv = this->mpParam->calculateCellSize(); + ParamVec cspv = mpParam->calculateCellSize(); if(LBMDIM==2) cspv[2] = 1.0; inmCellCnt = 1; double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0; @@ -1071,7 +1143,7 @@ bool LbmFsgrSolver::initializeSolverGrids() { //mStartSymm = false; #if ELBEEM_PLUGIN!=1 - if((LBMDIM==2)&&(this->mSizex<200)) { + if((LBMDIM==2)&&(mSizex<200)) { if(!checkSymmetry("init")) { errMsg("LbmFsgrSolver::initialize","Unsymmetric init..."); } else { @@ -1094,18 +1166,18 @@ bool LbmFsgrSolver::initializeSolverPostinit() { adaptGrid(lev); coarseRestrictFromFine(lev); } - this->markedClearList(); + markedClearList(); myTime_t fsgrtend = getTime(); - if(!this->mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<cDirNum; l++) { + for(int l=0; ldfVecX[l]!=0) area *= 0.5; - if(this->dfVecY[l]!=0) area *= 0.5; - if(this->dfVecZ[l]!=0) area *= 0.5; + if(dfVecX[l]!=0) area *= 0.5; + if(dfVecY[l]!=0) area *= 0.5; + if(dfVecZ[l]!=0) area *= 0.5; mFsgrCellArea[l] = area; } // l @@ -1123,6 +1195,8 @@ bool LbmFsgrSolver::initializeSolverPostinit() { stepMain(); // prepare once... + mpIso->setParticles(mpParticles, mPartDropMassSub); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<getSubdivs()<<", partsize="<mSmoothSurface<<","<mSmoothNormals<< /*","<mpIso->getSmoothSurface()<<","<mpIso->getSmoothNormals()<<") ",10); - debugOut("LbmFsgrSolver::initialize : Init done ... ",10); - this->mInitDone = 1; + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<getSmoothSurface()<<","<getSmoothNormals()<<") ",10); + debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10); + mInitDone = 1; #if LBM_INCLUDE_TESTSOLVERS==1 initCpdata(); @@ -1201,8 +1275,8 @@ bool LbmFsgrSolver::initializeSolverPostinit() { const LbmVec oldov=objvel; /*DEBUG*/ \ /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<mpGiObjects->size()); + int numobjs = (int)(mpGiObjects->size()); ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0; // 2d display as rectangles ntlVec3Gfx iniPos(0.0); if(LBMDIM==2) { dvec[2] = 1.0; - iniPos = (this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))-(dvec*0.0); + iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0); } else { - iniPos = (this->mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0); - //iniPos[2] = this->mvGeoStart[2] + dvec[2]*getForZMin1(); + iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0); } if( (int)mObjectMassMovnd.size() < numobjs) { @@ -1258,21 +1331,20 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0; int nbored; for(int OId=0; OIdmpGiObjects)[OId]; + ntlGeometryObject *obj = (*mpGiObjects)[OId]; bool skip = false; - if(obj->getGeoInitId() != this->mLbmInitId) skip=true; + if(obj->getGeoInitId() != mLbmInitId) skip=true; if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true; if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true; - debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<getName()<<" skip:"<getIsAnimated()<<" gid:"<getGeoInitId()<<" simgid:"<mLbmInitId, 10); if(skip) continue; + debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<getName()<<" skip:"<getIsAnimated()<<" gid:"<getGeoInitId()<<" simgid:"<getGeoInitType()&FGI_ALLBOUNDS) || (obj->getGeoInitType()&FGI_FLUID) && staticInit ) { otype = ntype = CFInvalid; switch(obj->getGeoInitType()) { - /* - case FGI_BNDPART: + /* case FGI_BNDPART: // old, use noslip for moving part/free objs case FGI_BNDFREE: if(!staticInit) { errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<getName() ); @@ -1311,8 +1383,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue; // copied from recalculateObjectSpeeds - mObjectSpeeds[OId] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ))); - debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 ); + mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ))); + debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 ); //vector tNormals; vector *pNormals = NULL; @@ -1324,7 +1396,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { // do two full update // TODO tNormals handling!? mMOIVerticesOld.clear(); - obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals, mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd); + obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals, mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd); monTrafo += mMOIVerticesOld.size(); obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false ); monTrafo += mMOIVertices.size(); @@ -1352,7 +1424,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { // FIXME? if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; } // get old type - CHECK FIXME , timestep could have changed - cause trouble? - ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - this->mpParam->getTimestep(),sourceTime); + ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime); if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; } } if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; } @@ -1360,6 +1432,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { LbmFloat massCheck = 0.; int massReinits=0; bool fillCells = (mObjectMassMovnd[OId]<=-1.); + LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime); + // first pass - set new obs. cells if(active) { @@ -1397,30 +1471,19 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { // compute fluid acceleration FORDF1 { if(rflagnb[l]&(CFFluid|CFInter)) { - const LbmFloat ux = this->dfDvecX[l]*objvel[0]; - const LbmFloat uy = this->dfDvecY[l]*objvel[1]; - const LbmFloat uz = this->dfDvecZ[l]*objvel[2]; + const LbmFloat ux = dfDvecX[l]*objvel[0]; + const LbmFloat uy = dfDvecY[l]*objvel[1]; + const LbmFloat uz = dfDvecZ[l]*objvel[2]; - /*LbmFloat *rhodstCell = RACPNT(level, i+dfVecX[l],j+dfVecY[l],k+dfVecZ[l],workSet); - LbmFloat rhodst = RAC(rhodstCell,dC); - for(int ll=1; ll<19; ll++) { rhodst += RAC(rhodstCell,ll); } - rhodst = 1.; // DEBUG TEST! */ - - LbmFloat factor = 2. * this->dfLength[l] * 3.0 * (ux+uy+uz); // - /* FSTEST - */ + LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // if(ntype&(CFBndFreeslip|CFBndPartslip)) { // missing, diag mass checks... - //if(l<=LBMDIM*2) massCheck += factor; - //if(l<=LBMDIM*2) factor *= 1.4142; - //if(l>LBMDIM*2) factor = 0.; - //else - //factor = 0.; // TODO, optimize factor *= 2.0; // TODO, optimize } else { factor *= 1.2; // TODO, optimize } + factor *= impactCorrFactor; // use manual tweaking channel RAC(dstCell,l) = factor; massCheck += factor; } else { @@ -1461,28 +1524,19 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) { } CellFlagType settype = CFInvalid; if(nbored&CFFluid) { - //if(nbored&(CFFluid|CFInter)) { settype = CFInter|CFNoInterpolSrc; rhomass = 1.5; if(!fillCells) rhomass = 0.; - //settype = CFFluid|CFNoInterpolSrc; rhomass=1.; // test - //rhomass = 1.01; // DEBUGT OBJVEL_CALC; if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; } - //settype=CFFluid; // rhomass = 1.; objvel = LbmVec(0.); // DEBUG TEST // new interpolate values LbmFloat avgrho = 0.0; LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0; interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz); - //objvel = LbmVec(avgux,avguy,avguz); - //avgrho=1.; initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) ); //errMsg("NMOCIT"," at "<mInitDone) { + if(mInitDone) { errMsg("initMov","dccd\n\nMassd test "<getName()<<" dccd massCheck="<0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<mLbmInitId <<") v"<initGeoTree(); - if(this->mAllfluid) { - this->freeGeoTree(); + initGeoTree(); + if(mAllfluid) { + freeGeoTree(); return true; } // make sure moving obstacles are inited correctly // preinit moving obj points - int numobjs = (int)(this->mpGiObjects->size()); + int numobjs = (int)(mpGiObjects->size()); for(int o=0; ompGiObjects)[o]; + ntlGeometryObject *obj = (*mpGiObjects)[o]; //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<getName()<<" type "<getGeoInitType()<<" anim"<getIsAnimated()<<" "<getVolumeInit() ,9); if( ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) || @@ -1677,20 +1729,20 @@ bool LbmFsgrSolver::initGeometryFlags() { } // max speed init - ntlVec3Gfx maxMovobjVelRw = this->getGeoMaxMovementVelocity( mSimulationTime, this->mpParam->getTimestep() ); - ntlVec3Gfx maxMovobjVel = vec2G( this->mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) ); - this->mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) ); - LbmFloat allowMax = this->mpParam->getTadapMaxSpeed(); // maximum allowed velocity + ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() ); + ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) ); + mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) ); + LbmFloat allowMax = mpParam->getTadapMaxSpeed(); // maximum allowed velocity debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<mpParam->setDesiredTimestep( newdt ); - this->mpParam->calculateAllMissingValues( mSimulationTime, this->mSilent ); - maxMovobjVel = vec2G( this->mpParam->calculateLattVelocityFromRw( vec2P(this->getGeoMaxMovementVelocity( - mSimulationTime, this->mpParam->getTimestep() )) )); + LbmFloat nextmax = mpParam->getSimulationMaxSpeed(); + LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr + debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5); + mpParam->setDesiredTimestep( newdt ); + mpParam->calculateAllMissingValues( mSimulationTime, mSilent ); + maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity( + mSimulationTime, mpParam->getTimestep() )) )); debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5); } recalculateObjectSpeeds(); @@ -1710,24 +1762,25 @@ bool LbmFsgrSolver::initGeometryFlags() { // 2d display as rectangles if(LBMDIM==2) { dvec[2] = 0.0; - iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))+(dvec*0.5); + iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5); //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } } } else { - iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5); + iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5); } // first init boundary conditions // invalid cells are set to empty afterwards + // TODO use floop macros!? for(int k= getForZMin1(); k< getForZMax1(level); ++k) { for(int j=1;jgeoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance); + const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance); if(inside) { - pObj = (*this->mpGiObjects)[OId]; + pObj = (*mpGiObjects)[OId]; switch( pObj->getGeoInitType() ){ case FGI_MBNDINFLOW: if(! pObj->getIsAnimated() ) { @@ -1795,7 +1848,7 @@ bool LbmFsgrSolver::initGeometryFlags() { ntype = CFInvalid; int inits = 0; GETPOS(i,j,k); - const bool inside = this->geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance); + const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance); if(inside) { ntype = CFFluid; } @@ -1838,13 +1891,15 @@ bool LbmFsgrSolver::initGeometryFlags() { } } - this->freeGeoTree(); + freeGeoTree(); myTime_t geotimeend = getTime(); debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; - if( TESTFLAG( RFLAG(mMaxRefine, ni, nj, nk, mLevel[mMaxRefine].setCurr), CFEmpty ) ) { + int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l]; + if( ( RFLAG(mMaxRefine, ni, nj, nk, mLevel[mMaxRefine].setCurr)& CFEmpty ) ) { LbmFloat arho=0., aux=0., auy=0., auz=0.; interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz); //errMsg("TINI"," "<cDirNum); l++) { - int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; + for(int l=0;(lcDirNum) ); + QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) ); QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass); } }}} @@ -2186,7 +2241,7 @@ bool LbmFsgrSolver::checkSymmetry(string idstring) LbmFloat maxdiv =0.0; if(erro) { errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv ); - //if(LBMDIM==2) this->mPanic = true; + //if(LBMDIM==2) mPanic = true; //return false; } else { errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv ); @@ -2205,17 +2260,17 @@ LbmFsgrSolver::interpolateCellValues( LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0; LbmFloat cellcnt = 0.0; - for(int nbl=1; nbl< this->cDfNum ; ++nbl) { + for(int nbl=1; nbl< cDfNum ; ++nbl) { if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue; if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){ //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) && //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) { cellcnt += 1.0; - for(int rl=0; rl< this->cDfNum ; ++rl) { + for(int rl=0; rl< cDfNum ; ++rl) { LbmFloat nbdf = QCELL_NB(level,ei,ej,ek, workSet,nbl, rl); - avgux += (this->dfDvecX[rl]*nbdf); - avguy += (this->dfDvecY[rl]*nbdf); - avguz += (this->dfDvecZ[rl]*nbdf); + avgux += (dfDvecX[rl]*nbdf); + avguy += (dfDvecY[rl]*nbdf); + avguz += (dfDvecZ[rl]*nbdf); avgrho += nbdf; } } @@ -2227,7 +2282,7 @@ LbmFsgrSolver::interpolateCellValues( avgux = avguy = avguz = 0.0; //TTT mNumProblems++; #if ELBEEM_PLUGIN!=1 - //this->mPanic=1; + //mPanic=1; // this can happen for worst case moving obj scenarios... errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<initId(); + + if(getenv("ELBEEM_RAWDEBUGDUMP")) { + debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2); + mDumpRawText = true; + } + + if(getenv("ELBEEM_BINDEBUGDUMP")) { + debMsgStd("LbmSolverInterface",DM_MSG,"Using env var ELBEEM_RAWDEBUGDUMP, mDumpRawText on",2); + mDumpRawBinary = true; + } } LbmSolverInterface::~LbmSolverInterface() @@ -79,8 +96,8 @@ void initGridSizes(int &sizex, int &sizey, int &sizez, int mMaxRefine, bool parallel) { // fix size inits to force cubic cells and mult4 level dimensions - const int debugGridsizeInit = 0; - if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<maxGridSize) maxGridSize = sizey; @@ -102,16 +119,19 @@ void initGridSizes(int &sizex, int &sizey, int &sizez, for(int i=0; i0.) { + memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) ); + } else { + // ignore 3 int slices... + memCnt += (double)( ( sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) ); + } if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"iso, mc:"<readString(name, "", source, target, needed); + string val = mpSifAttrs->readString(name, "", source, target, needed); if(!strcmp(val.c_str(),"")) { // no value given... return CFEmpty; @@ -212,12 +240,11 @@ CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultVal /*******************************************************************************/ /*! parse standard attributes */ void LbmSolverInterface::parseStdAttrList() { - if(!mpAttrs) { - errFatal("LbmSolverInterface::parseAttrList","mpAttrs pointer not initialized!",SIMWORLD_INITERROR); + if(!mpSifAttrs) { + errFatal("LbmSolverInterface::parseAttrList","mpSifAttrs pointer not initialized!",SIMWORLD_INITERROR); return; } // st currently unused - //mSurfaceTension = mpAttrs->readFloat("d_surfacetension", mSurfaceTension, "LbmSolverInterface", "mSurfaceTension", false); mBoundaryEast = readBoundaryFlagInt("boundary_east", mBoundaryEast, "LbmSolverInterface", "mBoundaryEast", false); mBoundaryWest = readBoundaryFlagInt("boundary_west", mBoundaryWest, "LbmSolverInterface", "mBoundaryWest", false); mBoundaryNorth = readBoundaryFlagInt("boundary_north", mBoundaryNorth,"LbmSolverInterface", "mBoundaryNorth", false); @@ -225,10 +252,10 @@ void LbmSolverInterface::parseStdAttrList() { mBoundaryTop = readBoundaryFlagInt("boundary_top", mBoundaryTop,"LbmSolverInterface", "mBoundaryTop", false); mBoundaryBottom= readBoundaryFlagInt("boundary_bottom", mBoundaryBottom,"LbmSolverInterface", "mBoundaryBottom", false); - mpAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo ); + mpSifAttrs->readMat4Gfx("domain_trafo" , (*mpSimTrafo), "ntlBlenderDumper","mpSimTrafo", false, mpSimTrafo ); LbmVec sizeVec(mSizex,mSizey,mSizez); - sizeVec = vec2L( mpAttrs->readVec3d("size", vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) ); + sizeVec = vec2L( mpSifAttrs->readVec3d("size", vec2P(sizeVec), "LbmSolverInterface", "sizeVec", false) ); mSizex = (int)sizeVec[0]; mSizey = (int)sizeVec[1]; mSizez = (int)sizeVec[2]; @@ -236,13 +263,13 @@ void LbmSolverInterface::parseStdAttrList() { // test solvers might not have mpPara, though if(mpParam) mpParam->setSize(mSizex, mSizey, mSizez ); - mInitDensityGradient = mpAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false); - mIsoValue = mpAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false ); + mInitDensityGradient = mpSifAttrs->readBool("initdensitygradient", mInitDensityGradient,"LbmSolverInterface", "mInitDensityGradient", false); + mIsoValue = mpSifAttrs->readFloat("isovalue", mIsoValue, "LbmOptSolver","mIsoValue", false ); - mDebugVelScale = mpAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false); - mNodeInfoString = mpAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false ); + mDebugVelScale = mpSifAttrs->readFloat("debugvelscale", mDebugVelScale,"LbmSolverInterface", "mDebugVelScale", false); + mNodeInfoString = mpSifAttrs->readString("nodeinfo", mNodeInfoString, "SimulationLbm","mNodeInfoString", false ); - mDumpVelocities = mpAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false ); + mDumpVelocities = mpSifAttrs->readBool("dump_velocities", mDumpVelocities, "SimulationLbm","mDumpVelocities", false ); if(getenv("ELBEEM_DUMPVELOCITIES")) { int get = atoi(getenv("ELBEEM_DUMPVELOCITIES")); if((get==0)||(get==1)) { @@ -252,14 +279,18 @@ void LbmSolverInterface::parseStdAttrList() { } // new test vars - mTForceStrength = 0.; // set from test solver mpAttrs->readFloat("tforcestrength", mTForceStrength,"LbmSolverInterface", "mTForceStrength", false); - mFarFieldSize = mpAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false); + // mTForceStrength set from test solver + mFarFieldSize = mpSifAttrs->readFloat("farfieldsize", mFarFieldSize,"LbmSolverInterface", "mFarFieldSize", false); // old compat - float sizeScale = mpAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false); + float sizeScale = mpSifAttrs->readFloat("test_scale", 0.,"LbmTestdata", "mSizeScale", false); if((mFarFieldSize<=0.)&&(sizeScale>0.)) { mFarFieldSize=sizeScale; errMsg("LbmTestdata","Warning - using mSizeScale..."); } - mCppfStage = mpAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false); - mPartGenProb = mpAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false); + mPartDropMassSub = mpSifAttrs->readFloat("part_dropmass", mPartDropMassSub,"LbmSolverInterface", "mPartDropMassSub", false); + + mCppfStage = mpSifAttrs->readInt("cppfstage", mCppfStage,"LbmSolverInterface", "mCppfStage", false); + mPartGenProb = mpSifAttrs->readFloat("partgenprob", mPartGenProb,"LbmFsgrSolver", "mPartGenProb", false); + mPartUsePhysModel = mpSifAttrs->readBool("partusephysmodel", mPartUsePhysModel,"LbmFsgrSolver", "mPartUsePhysModel", false); + mIsoSubdivs = mpSifAttrs->readInt("isosubdivs", mIsoSubdivs,"LbmFsgrSolver", "mIsoSubdivs", false); } @@ -278,9 +309,10 @@ void LbmSolverInterface::initGeoTree() { mGiObjSecondDist.resize( mpGiObjects->size() ); for(size_t i=0; isize(); i++) { if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true; - debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 ); + //debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 ); } debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<size() ,10); if(mpGiTree != NULL) delete mpGiTree; char treeFlag = (1<<(this->mLbmInitId+4)); @@ -307,15 +339,22 @@ int globGeoInitDebug = 0; int globGICPIProblems = 0; /*****************************************************************************/ /*! check for a certain flag type at position org */ -bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance) { +bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance, int shootDir) { // shift ve ctors to avoid rounding errors org += ntlVec3Gfx(0.0001); - ntlVec3Gfx dir = ntlVec3Gfx(1.0, 0.0, 0.0); OId = -1; + + // select shooting direction + ntlVec3Gfx dir = ntlVec3Gfx(1., 0., 0.); + if(shootDir==1) dir = ntlVec3Gfx(0., 1., 0.); + else if(shootDir==2) dir = ntlVec3Gfx(0., 0., 1.); + else if(shootDir==-2) dir = ntlVec3Gfx(0., 0., -1.); + else if(shootDir==-1) dir = ntlVec3Gfx(0., -1., 0.); + ntlRay ray(org, dir, 0, 1.0, mpGlob); - //int insCnt = 0; bool done = false; bool inside = false; + vector giObjFirstHistSide; giObjFirstHistSide.resize( mpGiObjects->size() ); for(size_t i=0; iintersectX(ray,distance,normal, triIns, flags, true); + if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true); + else mpGiTree->intersect(ray,distance,normal, triIns, flags, true); if(triIns) { ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance; LbmFloat orientation = dot(normal, dir); @@ -408,7 +448,8 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int ntlTriangle *triIns = NULL; distance = -1.0; ntlVec3Gfx normal(0.0); - mpGiTree->intersectX(ray,distance,normal, triIns, flags, true); + if(shootDir==0) mpGiTree->intersectX(ray,distance,normal, triIns, flags, true); + else mpGiTree->intersect(ray,distance,normal, triIns, flags, true); if(triIns) { // check outside intersect LbmFloat orientation = dot(normal, dir); diff --git a/intern/elbeem/intern/solver_interface.h b/intern/elbeem/intern/solver_interface.h index b95ce0358cb..e3bc06a284c 100644 --- a/intern/elbeem/intern/solver_interface.h +++ b/intern/elbeem/intern/solver_interface.h @@ -3,7 +3,7 @@ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method * All code distributed as part of El'Beem is covered by the version 2 of the * GNU General Public License. See the file COPYING for details. - * Copyright 2003-2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Header for Combined 2D/3D Lattice Boltzmann Interface Class * @@ -33,6 +33,7 @@ #include "attributes.h" #include "isosurface.h" class ParticleTracer; +class ParticleObject; // use which fp-precision for LBM? 1=float, 2=double #ifdef PRECISION_LBM_SINGLE @@ -135,12 +136,13 @@ typedef int BubbleId; #define CFInvalid (CellFlagType)(1<<31) // use 32bit flag types -#ifdef __x86_64__ - typedef int cfINT32; -#else - typedef long cfINT32; -#endif // defined (_IA64) -#define CellFlagType cfINT32 +//#ifdef __x86_64__ + //typedef int cfINT32; +//#else + //typedef long cfINT32; +//#endif // defined (_IA64) +//#define CellFlagType cfINT32 +#define CellFlagType int #define CellFlagTypeSize 4 @@ -261,7 +263,7 @@ class LbmSolverInterface virtual int initParticles() = 0; virtual void advanceParticles() = 0; /*! get surface object (NULL if no surface) */ - ntlGeometryObject* getSurfaceGeoObj() { return mpIso; } + IsoSurface* getSurfaceGeoObj() { return mpIso; } /*! debug object display */ virtual vector getDebugObjects() { vector empty(0); return empty; } @@ -276,7 +278,7 @@ class LbmSolverInterface /*! destroy tree etc. when geometry init done */ void freeGeoTree(); /*! check for a certain flag type at position org (needed for e.g. quadtree refinement) */ - bool geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance); + bool geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance,int shootDir=0); bool geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse); /*! set render globals, for scene/tree access */ @@ -306,9 +308,12 @@ class LbmSolverInterface bool getAllfluid() { return mAllfluid; } /*! set attr list pointer */ - void setAttrList(AttributeList *set) { mpAttrs = set; } + void setAttrList(AttributeList *set) { mpSifAttrs = set; } /*! Returns the attribute list pointer */ - inline AttributeList *getAttributeList() { return mpAttrs; } + inline AttributeList *getAttributeList() { return mpSifAttrs; } + /*! set sws attr list pointer */ + void setSwsAttrList(AttributeList *set) { mpSifSwsAttrs = set; } + inline AttributeList *getSwsAttributeList() { return mpSifSwsAttrs; } /*! set parametrizer pointer */ inline void setParametrizer(Parametrizer *set) { mpParam = set; } @@ -353,6 +358,8 @@ class LbmSolverInterface //! set amount of surface/normal smoothing inline void setSmoothing(float setss,float setns){ mSmoothSurface=setss; mSmoothNormals=setns; } + //! set amount of iso subdivisions + inline void setIsoSubdivs(int s){ mIsoSubdivs=s; } //! set desired refinement inline void setPreviewSize(int set){ mOutputSurfacePreview = set; } //! set desired refinement @@ -378,6 +385,13 @@ class LbmSolverInterface //! set/get cp stage inline void setCpStage(int set) { mCppfStage = set; } inline int getCpStage() const { return mCppfStage; } + //! set/get dump modes + inline void setDumpRawText(bool set) { mDumpRawText = set; } + inline bool getDumpRawText() const { return mDumpRawText; } + inline void setDumpRawBinary(bool set) { mDumpRawBinary = set; } + inline bool getDumpRawBinary() const { return mDumpRawBinary; } + inline void setDumpRawBinaryZip(bool set) { mDumpRawBinaryZip = set; } + inline bool getDumpRawBinaryZip() const { return mDumpRawBinaryZip; } // cell iterator interface @@ -468,8 +482,9 @@ class LbmSolverInterface /*! init density gradient? */ bool mInitDensityGradient; - /*! pointer to the attribute list */ - AttributeList *mpAttrs; + /*! pointer to the attribute list, only pointer to obj attrs */ + AttributeList *mpSifAttrs; + AttributeList *mpSifSwsAttrs; /*! get parameters from this parametrize in finishInit */ Parametrizer *mpParam; @@ -534,9 +549,11 @@ class LbmSolverInterface int mOutputSurfacePreview; LbmFloat mPreviewFactor; - /* enable surface and normals smoothing? */ + /*! enable surface and normals smoothing? */ float mSmoothSurface; float mSmoothNormals; + /*! isosurface subdivisions */ + int mIsoSubdivs; //! particle generation probability LbmFloat mPartGenProb; @@ -553,14 +570,22 @@ class LbmSolverInterface //! part slip value for domain LbmFloat mDomainPartSlipValue; + // size of far field area + LbmFloat mFarFieldSize; + // amount of drop mass to subtract + LbmFloat mPartDropMassSub; + // use physical drop model for particles? + bool mPartUsePhysModel; + //! test vars // strength of applied force LbmFloat mTForceStrength; - // - LbmFloat mFarFieldSize; - // int mCppfStage; + //! dumping modes + bool mDumpRawText; + bool mDumpRawBinary; + bool mDumpRawBinaryZip; }; @@ -569,7 +594,7 @@ void initGridSizes(int &mSizex, int &mSizey, int &mSizez, ntlVec3Gfx &mvGeoStart, ntlVec3Gfx &mvGeoEnd, int mMaxRefine, bool parallel); void calculateMemreqEstimate(int resx,int resy,int resz, int refine, - double *reqret, string *reqstr); + float farfieldsize, double *reqret, string *reqstr); //! helper function to convert flag to string (for debuggin) string convertCellFlagType2String( CellFlagType flag ); diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp index 964192420e4..bced75ab444 100644 --- a/intern/elbeem/intern/solver_main.cpp +++ b/intern/elbeem/intern/solver_main.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004,2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Standard LBM Factory implementation * @@ -20,65 +20,30 @@ double globdfcnt; double globdfavg[19]; double globdfmax[19]; double globdfmin[19]; +extern int glob_mpindex,glob_mpnum; +extern bool globOutstrForce; +// simulation object interface void LbmFsgrSolver::step() { - initLevelOmegas(); - stepMain(); - - // intlbm test - if(0) { - if(this->mStepCnt<5) { - // init - globdfcnt=0.; - FORDF0{ - globdfavg[l] = 0.; - globdfmax[l] = -1000.; //this->dfEquil[l]; - globdfmin[l] = 1000.; // this->dfEquil[l]; - } - } else { - - int workSet = mLevel[mMaxRefine].setCurr; - for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { - for(int j=0;jglobdfmax[l]) globdfmax[l] = df; - //if(df<=0.01) { errMsg("GLOBDFERR"," at "<mStepCnt%10==0) { - FORDF0{ errMsg("GLOBDF","l="<markedClearList(); // DMC clearMarkedCellsList +// lbm main step +void messageOutputForce(string from); +void LbmFsgrSolver::stepMain() { + myTime_t timestart = getTime(); + + initLevelOmegas(); + markedClearList(); // DMC clearMarkedCellsList // safety check, counter reset - this->mNumUsedCells = 0; + mNumUsedCells = 0; mNumInterdCells = 0; mNumInvIfCells = 0; - //debugOutNnl("LbmFsgrSolver::step : "<mStepCnt, 10); - if(!this->mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<mStepCnt<<" t:"<mStepCnt<<" "); - myTime_t timestart = getTime(); + //debugOutNnl("LbmFsgrSolver::step : "<mFixMass; // reset here for next step + mCurrentMass = mFixMass; // reset here for next step mCurrentVolume = 0.0; //change to single step advance! int levsteps = 0; - int dsbits = this->mStepCnt ^ (this->mStepCnt-1); - //errMsg("S"," step:"<mStepCnt<<" s-1:"<<(this->mStepCnt-1)<<" xf:"<mStepCnt&(1<mStepCnt++; + mStepCnt++; // some dbugging output follows // calculate MLSUPS myTime_t timeend = getTime(); - this->mNumUsedCells += mNumInterdCells; // count both types for MLSUPS - mAvgNumUsedCells += this->mNumUsedCells; - this->mMLSUPS = (this->mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000); - if(this->mMLSUPS>10000){ this->mMLSUPS = -1; } - else { mAvgMLSUPS += this->mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups + mNumUsedCells += mNumInterdCells; // count both types for MLSUPS + mAvgNumUsedCells += mNumUsedCells; + mMLSUPS = ((double)mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000.0); + if(mMLSUPS>10000){ mMLSUPS = -1; } + //else { mAvgMLSUPS += mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000); if(totMLSUPS>10000) totMLSUPS = -1; mNumInvIfTotal += mNumInvIfCells; // debug // do some formatting - if(!this->mSilent){ - string sepStr(""); // DEBUG - int avgcls = (int)(mAvgNumUsedCells/(LONGINT)this->mStepCnt); - debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<mStepCnt<<" t:"<mMLSUPS<< - " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<< - " totcls:"<<(this->mNumUsedCells)<< sepStr<< - " avgcls:"<< avgcls<< sepStr<< - " intd:"<mNumFilledCells<<", emptied:"<mNumEmptiedCells<< sepStr<< - " mMxv:"<mName<<"' " , 10); + if(!mSilent){ + int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt); + debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<mStepCnt==1) { - mMinNoCells = mMaxNoCells = this->mNumUsedCells; + if(mStepCnt==1) { + mMinNoCells = mMaxNoCells = mNumUsedCells; } else { - if(this->mNumUsedCells>mMaxNoCells) mMaxNoCells = this->mNumUsedCells; - if(this->mNumUsedCellsmNumUsedCells; + if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells; + if(mNumUsedCellsmStepCnt<<" levstep:"<mInitDone)) { handleTestdata(); } - testXchng(); + if((mUseTestdata)&&(mInitDone)) { handleTestdata(); } + mrExchange(); #endif // advance positions with current grid @@ -258,18 +217,28 @@ void LbmFsgrSolver::stepMain() #endif // WIN32 // output total step time - timeend = getTime(); + myTime_t timeend2 = getTime(); double mdelta = (lastMass-mCurrentMass); if(ABS(mdelta)<1e-12) mdelta=0.; - debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<mStepCnt - <<": dccd="<< mCurrentMass - <<",d"<10000){ effMLSUPS = -1; } + else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups + } + + debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<mpParam->setFluidVolumeHeight(mFVHeight); + mpParam->setFluidVolumeHeight(mFVHeight); } // advance time before timestep change - mSimulationTime += this->mpParam->getTimestep(); + mSimulationTime += mpParam->getTimestep(); // time adaptivity - this->mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) ); + mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) ); //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG - if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); } + if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); } // DEBUG VEL CHECK if(0) { @@ -370,10 +339,10 @@ void LbmFsgrSolver::fineAdvance() #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5) #if LBMDIM==3 - // normal -#define SLOWDOWNREGION (this->mSizez/4) +// normal +#define SLOWDOWNREGION (mSizez/4) #else // LBMDIM==2 - // off +// off #define SLOWDOWNREGION 10 #endif // LBMDIM==2 #define P_LCSMQO 0.01 @@ -397,7 +366,7 @@ LbmFsgrSolver::mainLoop(int lev) # if LBM_INCLUDE_TESTSOLVERS==1 // 3d region off... quit - if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; } + if((mUseTestdata)&&(mpTest->mFarfMode>0)) { return; } # endif // ELBEEM_PLUGIN!=1 // main loop region @@ -448,16 +417,16 @@ LbmFsgrSolver::mainLoop(int lev) errMsg("LbmFsgrSolver::mainLoop","Err flagp "<mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) { + if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) { FORDF0 { RAC(tcel,l) = RAC(ccel,l); } } else { interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false); @@ -731,27 +700,14 @@ LbmFsgrSolver::mainLoop(int lev) } // l // normal interface, no if empty/fluid - LbmFloat nv1,nv2; - LbmFloat nx,ny,nz; + // computenormal + LbmFloat surfaceNormal[3]; + computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal); - if(nbflag[dE] &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0; - if(nbflag[dW] &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0; - nx = 0.5* (nv2-nv1); - if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0; - if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0; - ny = 0.5* (nv2-nv1); -# if LBMDIM==3 - if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0; - if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0; - nz = 0.5* (nv2-nv1); -# else // LBMDIM==3 - nz = 0.0; -# endif // LBMDIM==3 - - if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) { + if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) { // normal ok and usable... FORDF1 { - if( (this->dfDvecX[l]*nx + this->dfDvecY[l]*ny + this->dfDvecZ[l]*nz) // dot Dvec,norml + if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2]) // dot Dvec,norml > LBM_EPSILON) { recons[l] = 2; numRecons++; @@ -781,7 +737,8 @@ LbmFsgrSolver::mainLoop(int lev) - MYDF( l ); } } - usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on + ux=oldUx, uy=oldUy, uz=oldUz; // no local vars, only for usqr + usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on # else // OPT3D==0 oldRho = + RAC(ccel,dC) + RAC(ccel,dN ) + RAC(ccel,dS ) + RAC(ccel,dE ) @@ -864,11 +821,15 @@ LbmFsgrSolver::mainLoop(int lev) } FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); } } else {// NEWSURFT */ + if(usqr>0.3*0.3) { + // force reset! , warning causes distortions... + FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); } + } else { // normal collide // mass streaming done... do normal collide LbmVec grav = mLevel[lev].gravity*mass; DEFAULT_COLLIDEG(grav); - PERFORM_USQRMAXCHECK; + PERFORM_USQRMAXCHECK; } // rho init from default collide necessary for fill/empty check below } // test } @@ -876,87 +837,111 @@ LbmFsgrSolver::mainLoop(int lev) // testing..., particle generation // also check oldFlag for CFNoNbFluid, test // for inflow no pargen test // NOBUBBB! - if((this->mInitDone) //&&(mUseTestdata) - //&& (!((oldFlag|newFlag)&CFNoNbEmpty)) - && (!((oldFlag|newFlag)&CFNoDelete)) - && (this->mPartGenProb>0.0)) { + if((mInitDone) + // dont allow new if cells, or submerged ones + && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) + // dont try to subtract from empty cells + && (mass>0.) && (mPartGenProb>0.0)) { bool doAdd = true; bool bndOk=true; - if( (ithis->mSizex-cutMin)|| - (jthis->mSizey-cutMin)|| - (kthis->mSizez-cutMin) ) { bndOk=false; } + if( (imSizex-cutMin)|| + (jmSizey-cutMin)|| + (kmSizez-cutMin) ) { bndOk=false; } if(!bndOk) doAdd=false; - LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep); - LbmFloat rux = (ux * realWorldFac); - LbmFloat ruy = (uy * realWorldFac); - LbmFloat ruz = (uz * realWorldFac); - LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz)); - // WHMOD - LbmFloat prob = (rand()/(RAND_MAX+1.0)); - LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl; + LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333; - // reduce probability in outer region? - const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst; - const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst; - LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord); - LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord); - if(pifac<0.) pifac=0.; - if(pjfac<0.) pjfac=0.; + // physical drop model + if(mPartUsePhysModel) { + LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep); + LbmFloat rux = (ux * realWorldFac); + LbmFloat ruy = (uy * realWorldFac); + LbmFloat ruz = (uz * realWorldFac); + LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz)); + basethresh *= rl; - //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { - if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { - // add - } else { - doAdd = false; // dont... - } + // reduce probability in outer region? + const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst; + const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst; + LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord); + LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord); + if(pifac<0.) pifac=0.; + if(pjfac<0.) pjfac=0.; - // "wind" disturbance - // use realworld relative velocity here instead? - if( (doAdd && - ((rl>RWVEL_WINDTHRESH) && (lcsmqothis->mSizez-SLOWDOWNREGION) ) { - LbmFloat nuz = uz; - if(k>this->mSizez-SLOWDOWNREGION) { - // special case - LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) ); - zfac /= (LbmFloat)(SLOWDOWNREGION); - nuz += (1.0) * zfac; // check max speed? OFF? - //errMsg("TOPT"," at "<0.0095) && (rl>RWVEL_THRESH) ) { + if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { + // add } else { - // normal probability - //? LbmFloat fac = P_LCSMQO-lcsmqo; - //? jdf *= fac; + doAdd = false; // dont... } - FORDF1 { - const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0)); - // TODO use wind velocity? - if(jdf>0.025) { - const LbmFloat add = this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf; - RAC(tcel,l) += add; } - } - //errMsg("TOPDOWNCORR"," jdf:"<RWVEL_WINDTHRESH) && (lcsmqomSizez-SLOWDOWNREGION) ) { + LbmFloat nuz = uz; + if(k>mSizez-SLOWDOWNREGION) { + // special case + LbmFloat zfac = (LbmFloat)( k-(mSizez-SLOWDOWNREGION) ); + zfac /= (LbmFloat)(SLOWDOWNREGION); + nuz += (1.0) * zfac; // check max speed? OFF? + //errMsg("TOPT"," at "<0.025) { + const LbmFloat add = this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf; + RAC(tcel,l) += add; } + } + //errMsg("TOPDOWNCORR"," jdf:"<0.0095)) { // add + if((prob0.012)) { // add + } else { doAdd = false; }// dont... + } + + + // remove noise 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(this->mStepCnt>700) errMsg("DFJITT"," at "<getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel() ); + // subtract drop mass + mass -= dropmass; + // init new particle + { + ParticleObject np( ntlVec3Gfx(srci,srcj,srck) ); + np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]); + np.setStatus(PART_IN); + np.setType(type); + //if((s%3)==2) np.setType(PART_FLOAT); + np.setSize(newPartsize); + //errMsg("NEWPART"," at "<mNumFilledCells++; // DEBUG + //mNumFilledCells++; // DEBUG calcCellsFilled++; } else if(ifemptied) { @@ -1055,7 +1046,7 @@ LbmFsgrSolver::mainLoop(int lev) if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; // NEWSURFT emptyp.x = i; emptyp.y = j; emptyp.z = k; LIST_EMPTY(emptyp); - //this->mNumEmptiedCells++; // DEBUG + //mNumEmptiedCells++; // DEBUG calcCellsEmptied++; } // dont cutoff values -> better cell conversions @@ -1092,13 +1083,6 @@ LbmFsgrSolver::mainLoop(int lev) // interface cell handling done... -//#if PARALLEL==1 -//#include "paraloopend.h" - //GRID_REGION_END(); -//#else // PARALLEL==1 - //GRID_LOOPREG_END(); - //GRID_REGION_END(); -//#endif // PARALLEL==1 #if PARALLEL!=1 GRID_LOOPREG_END(); #else // PARALLEL==1 @@ -1108,9 +1092,9 @@ LbmFsgrSolver::mainLoop(int lev) // write vars from computations to class mLevel[lev].lmass = calcCurrentMass; mLevel[lev].lvolume = calcCurrentVolume; - this->mNumFilledCells = calcCellsFilled; - this->mNumEmptiedCells = calcCellsEmptied; - this->mNumUsedCells = calcNumUsedCells; + mNumFilledCells = calcCellsFilled; + mNumEmptiedCells = calcCellsEmptied; + mNumUsedCells = calcNumUsedCells; } @@ -1122,34 +1106,32 @@ LbmFsgrSolver::preinitGrids() const bool doReduce = false; const int gridLoopBound=0; - // touch both grids + // preinit both grids for(int s=0; s<2; s++) { - GRID_REGION_INIT(); + GRID_REGION_INIT(); #if PARALLEL==1 #include "paraloopstart.h" #endif // PARALLEL==1 - GRID_REGION_START(); - GRID_LOOP_START(); - FORDF0{ RAC(ccel,l) = 0.; } - *pFlagSrc =0; - *pFlagDst =0; - //errMsg("l1"," at "<=mSizex-1) { nx = -1.; nonorm = true; } + if(j<=0) { ny = 1.; nonorm = true; } + if(j>=mSizey-1) { ny = -1.; nonorm = true; } +# if LBMDIM==3 + if(k<=0) { nz = 1.; nonorm = true; } + if(k>=mSizez-1) { nz = -1.; nonorm = true; } +# endif // LBMDIM==3 + if(!nonorm) { + // in domain, revert to helper, use setCurr&mMaxRefine + LbmVec bnormal; + CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr); + LbmFloat *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr); + computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]); + // TODO check if there is a normal near here? + // use wider range otherwise... + snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0]; + return; + } + snret[0]=nx; snret[1]=ny; snret[2]=nz; +} + void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) { #if FSGR_STRICT_DEBUG==10 // dangerous, this can change the simulation... @@ -1494,9 +1567,9 @@ void LbmFsgrSolver::reinitFlags( int workSet ) { massChange = 0.0; } else { // Problem! no interface neighbors - this->mFixMass += massChange; + mFixMass += massChange; //TTT mNumProblems++; - //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<mFixMass); + //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<mFixMass += massChange; + mFixMass += massChange; //TTT mNumProblems++; - //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<mFixMass); + //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<mFixMass * newIfFac); + QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac); int nbored = 0; FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); } @@ -1602,15 +1675,16 @@ void LbmFsgrSolver::reinitFlags( int workSet ) { int i=iter->x, j=iter->y, k=iter->z; if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; } - LbmFloat nrho = 0.0; - FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); } - QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho; - QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT; + initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) { + //LbmFloat nrho = 0.0; + //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); } + //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho; + //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT; } if(mListNewInter.size()>0){ - //errMsg("FixMassDisted"," fm:"<mFixMass<<" nif:"<mFixMass = 0.0; + //errMsg("FixMassDisted"," fm:"<dfVecX[this->dfInv[l]], nj=j+this->dfVecY[this->dfInv[l]], nk=k+this->dfVecZ[this->dfInv[l]]; \ + errMsg("STREAMCHECK"," Invalid nbflag, streamed DF l"<dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \ } \ } \ @@ -1198,5 +1207,7 @@ inline void LbmFsgrSolver::collideArrays( muy = uy; muz = uz; outrho = rho; + + lev=i=j=k; // debug, remove warnings }; diff --git a/intern/elbeem/intern/solver_util.cpp b/intern/elbeem/intern/solver_util.cpp index 6b7f6f2ed9c..d11ae6c6b40 100644 --- a/intern/elbeem/intern/solver_util.cpp +++ b/intern/elbeem/intern/solver_util.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004,2005 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Standard LBM Factory implementation * @@ -15,6 +15,7 @@ #include "ntl_world.h" #include "simulation_object.h" +#include /****************************************************************************** * helper functions @@ -23,6 +24,10 @@ // try to enhance surface? #define SURFACE_ENH 2 +extern bool glob_mpactive; +extern bool glob_mpnum; +extern bool glob_mpindex; + //! for raytracing void LbmFsgrSolver::prepareVisualization( void ) { int lev = mMaxRefine; @@ -36,10 +41,8 @@ void LbmFsgrSolver::prepareVisualization( void ) { mainGravLen = thisGravLen; mainGravDir = l; } - //errMsg("GRAVT","l"<mpIso->lbmGetData(i,j,ZKOFF)=0.0; + *mpIso->lbmGetData(i,j,ZKOFF)=0.0; } #else // LBMDIM==2 // 3d, use normal bounds @@ -59,21 +62,17 @@ void LbmFsgrSolver::prepareVisualization( void ) { for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) for(int j=0;jmpIso->lbmGetData(i,j,ZKOFF)=0.0; + *mpIso->lbmGetData(i,j,ZKOFF)=0.0; } #endif // LBMDIM==2 // MPT, ignore - //if( strstr( this->getName().c_str(), "mpfluid2" ) != NULL) { - //errMsg("MPT TESTX","skipping mpfluid2"); - //return; - //} - if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) { + if((glob_mpactive) && (glob_mpnum>1) && (glob_mpindex==0)) { mpIso->resetAll(0.); } - LbmFloat minval = this->mIsoValue*1.05; // / mIsoWeight[13]; + LbmFloat minval = mIsoValue*1.05; // / mIsoWeight[13]; // add up... float val = 0.0; for(int k= getForZMin1(); k< getForZMax1(lev); ++k) @@ -117,10 +116,10 @@ void LbmFsgrSolver::prepareVisualization( void ) { if(cflag&CFEmpty) { // special empty treatment if((noslipbnd)&&(intercnt>6)) { - *this->mpIso->lbmGetData(i,j,ZKOFF) += minval; + *mpIso->lbmGetData(i,j,ZKOFF) += minval; } else if((noslipbnd)&&(intercnt>0)) { // necessary? - *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.9; + *mpIso->lbmGetData(i,j,ZKOFF) += mIsoValue*0.9; } else { // nothing to do... } @@ -135,7 +134,7 @@ void LbmFsgrSolver::prepareVisualization( void ) { if(noslipbnd) { if(valmpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); + *mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); } #endif // SURFACE_ENH>0 @@ -143,43 +142,43 @@ void LbmFsgrSolver::prepareVisualization( void ) { continue; } - *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); - *this->mpIso->lbmGetData( i , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); - *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); + *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); + *mpIso->lbmGetData( i , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); + *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); - *this->mpIso->lbmGetData( i-1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); - *this->mpIso->lbmGetData( i , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); - *this->mpIso->lbmGetData( i+1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); + *mpIso->lbmGetData( i-1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); + *mpIso->lbmGetData( i , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); + *mpIso->lbmGetData( i+1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); - *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); - *this->mpIso->lbmGetData( i , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); - *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); + *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); + *mpIso->lbmGetData( i , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); + *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); - *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[9] ); - *this->mpIso->lbmGetData( i , j-1 ,ZKOFF ) += ( val * mIsoWeight[10] ); - *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[11] ); + *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[9] ); + *mpIso->lbmGetData( i , j-1 ,ZKOFF ) += ( val * mIsoWeight[10] ); + *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[11] ); - *this->mpIso->lbmGetData( i-1 , j ,ZKOFF ) += ( val * mIsoWeight[12] ); - *this->mpIso->lbmGetData( i , j ,ZKOFF ) += ( val * mIsoWeight[13] ); - *this->mpIso->lbmGetData( i+1 , j ,ZKOFF ) += ( val * mIsoWeight[14] ); + *mpIso->lbmGetData( i-1 , j ,ZKOFF ) += ( val * mIsoWeight[12] ); + *mpIso->lbmGetData( i , j ,ZKOFF ) += ( val * mIsoWeight[13] ); + *mpIso->lbmGetData( i+1 , j ,ZKOFF ) += ( val * mIsoWeight[14] ); - *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[15] ); - *this->mpIso->lbmGetData( i , j+1 ,ZKOFF ) += ( val * mIsoWeight[16] ); - *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[17] ); + *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[15] ); + *mpIso->lbmGetData( i , j+1 ,ZKOFF ) += ( val * mIsoWeight[16] ); + *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[17] ); - *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); - *this->mpIso->lbmGetData( i , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); - *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); + *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); + *mpIso->lbmGetData( i , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); + *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); - *this->mpIso->lbmGetData( i-1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); - *this->mpIso->lbmGetData( i , j ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); - *this->mpIso->lbmGetData( i+1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); + *mpIso->lbmGetData( i-1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); + *mpIso->lbmGetData( i , j ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); + *mpIso->lbmGetData( i+1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); - *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); - *this->mpIso->lbmGetData( i , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); - *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); + *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); + *mpIso->lbmGetData( i , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); + *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); } // TEST!? @@ -194,9 +193,9 @@ void LbmFsgrSolver::prepareVisualization( void ) { CellFlagType nbored=0; LbmFloat avgfill=0.,avgfcnt=0.; FORDF1 { - const int ni = i+this->dfVecX[l]; - const int nj = j+this->dfVecY[l]; - const int nk = ZKOFF+this->dfVecZ[l]; + const int ni = i+dfVecX[l]; + const int nj = j+dfVecY[l]; + const int nk = ZKOFF+dfVecZ[l]; const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet); nbored |= nbflag; @@ -216,10 +215,10 @@ void LbmFsgrSolver::prepareVisualization( void ) { if(nbored&CFInter) { if(avgfcnt>0.) avgfill/=avgfcnt; - *this->mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue; + *mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue; } else if(nbored&CFFluid) { - *this->mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue; + *mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue; } } @@ -232,16 +231,16 @@ void LbmFsgrSolver::prepareVisualization( void ) { for(int j=1;jmpIso->lbmGetData(i,j,ZKOFF)==0.)) { + if( (cflag&(CFBnd)) && (*mpIso->lbmGetData(i,j,ZKOFF)==0.)) { int bndnbcnt=0; FORDF1 { - const int ni = i+this->dfVecX[l]; - const int nj = j+this->dfVecY[l]; - const int nk = ZKOFF+this->dfVecZ[l]; + const int ni = i+dfVecX[l]; + const int nj = j+dfVecY[l]; + const int nk = ZKOFF+dfVecZ[l]; const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet); if(nbflag&CFBnd) bndnbcnt++; } - if(bndnbcnt>0) *this->mpIso->lbmGetData(i,j,ZKOFF)=this->mIsoValue*0.95; + if(bndnbcnt>0) *mpIso->lbmGetData(i,j,ZKOFF)=mIsoValue*0.95; } } } @@ -250,70 +249,70 @@ void LbmFsgrSolver::prepareVisualization( void ) { if(mFsSurfGenSetting&fssgNoNorth) for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) for(int j=0;jmpIso->lbmGetData(0, j,ZKOFF) = *this->mpIso->lbmGetData(1, j,ZKOFF); + *mpIso->lbmGetData(0, j,ZKOFF) = *mpIso->lbmGetData(1, j,ZKOFF); } if(mFsSurfGenSetting&fssgNoEast) for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) for(int i=0;impIso->lbmGetData(i,0, ZKOFF) = *this->mpIso->lbmGetData(i,1, ZKOFF); + *mpIso->lbmGetData(i,0, ZKOFF) = *mpIso->lbmGetData(i,1, ZKOFF); } if(mFsSurfGenSetting&fssgNoSouth) for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) for(int j=0;jmpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *this->mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF); + *mpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF); } if(mFsSurfGenSetting&fssgNoWest) for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) for(int i=0;impIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF); + *mpIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF); } if(LBMDIM>2) { if(mFsSurfGenSetting&fssgNoBottom) for(int j=0;jmpIso->lbmGetData(i,j,0 ) = *this->mpIso->lbmGetData(i,j,1 ); + *mpIso->lbmGetData(i,j,0 ) = *mpIso->lbmGetData(i,j,1 ); } if(mFsSurfGenSetting&fssgNoTop) for(int j=0;jmpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2); + *mpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2); } } #endif // SURFACE_ENH>=2 // update preview, remove 2d? - if((this->mOutputSurfacePreview)&&(LBMDIM==3)) { - int pvsx = (int)(this->mPreviewFactor*this->mSizex); - int pvsy = (int)(this->mPreviewFactor*this->mSizey); - int pvsz = (int)(this->mPreviewFactor*this->mSizez); - //float scale = (float)this->mSizex / previewSize; - LbmFloat scalex = (LbmFloat)this->mSizex/(LbmFloat)pvsx; - LbmFloat scaley = (LbmFloat)this->mSizey/(LbmFloat)pvsy; - LbmFloat scalez = (LbmFloat)this->mSizez/(LbmFloat)pvsz; + if((mOutputSurfacePreview)&&(LBMDIM==3)) { + int pvsx = (int)(mPreviewFactor*mSizex); + int pvsy = (int)(mPreviewFactor*mSizey); + int pvsz = (int)(mPreviewFactor*mSizez); + //float scale = (float)mSizex / previewSize; + LbmFloat scalex = (LbmFloat)mSizex/(LbmFloat)pvsx; + LbmFloat scaley = (LbmFloat)mSizey/(LbmFloat)pvsy; + LbmFloat scalez = (LbmFloat)mSizez/(LbmFloat)pvsz; for(int k= 0; k< (pvsz-1); ++k) for(int j=0;j< pvsy;j++) for(int i=0;i< pvsx;i++) { - *mpPreviewSurface->lbmGetData(i,j,k) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) ); + *mpPreviewSurface->lbmGetData(i,j,k) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) ); } // set borders again... for(int k= 0; k< (pvsz-1); ++k) { for(int j=0;j< pvsy;j++) { - *mpPreviewSurface->lbmGetData(0,j,k) = *this->mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) ); - *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *this->mpIso->lbmGetData( this->mSizex-1, (int)(j*scaley), (int)(k*scalez) ); + *mpPreviewSurface->lbmGetData(0,j,k) = *mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) ); + *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *mpIso->lbmGetData( mSizex-1, (int)(j*scaley), (int)(k*scalez) ); } for(int i=0;i< pvsx;i++) { - *mpPreviewSurface->lbmGetData(i,0,k) = *this->mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) ); - *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *this->mpIso->lbmGetData( (int)(i*scalex), this->mSizey-1, (int)(k*scalez) ); + *mpPreviewSurface->lbmGetData(i,0,k) = *mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) ); + *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *mpIso->lbmGetData( (int)(i*scalex), mSizey-1, (int)(k*scalez) ); } } for(int j=0;jlbmGetData(i,j,0) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0); - *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , this->mSizez-1); + *mpPreviewSurface->lbmGetData(i,j,0) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0); + *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , mSizez-1); } - if(mFarFieldSize>=2.) { + if(mFarFieldSize>=1.2) { // also remove preview border for(int k= 0; k< (pvsz-1); ++k) { for(int j=0;j< pvsy;j++) { @@ -349,45 +348,9 @@ void LbmFsgrSolver::prepareVisualization( void ) { } // MPT -#if LBM_INCLUDE_TESTSOLVERS==1 - if( ( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) && (mInitDone) ) { - errMsg("MPT TESTX","special mpfluid1"); - vector *sims = mpGlob->getSims(); - for(int simi=0; simi<(int)(sims->size()); simi++) { - SimulationObject *sim = (*sims)[simi]; - if( strstr( sim->getName().c_str(), "mpfluid2" ) != NULL) { - LbmFsgrSolver *fsgr = (LbmFsgrSolver *)(sim->getSolver()); - int msyLev = mMaxRefine; - int simLev = fsgr->mMaxRefine; - - IsoSurface* simIso = fsgr->mpIso; - int idst = mLevel[msyLev].lSizex-2 -1; // start at boundary - int isrc = 0 +2; - if((simIso)&&(fsgr->mInitDone)) { - fsgr->prepareVisualization(); - - for(int k=0;k<=mLevel[simLev].lSizez-1;++k) { - for(int j=0;j<=mLevel[simLev].lSizey-1;++j) { - for(int i=isrc; ilbmGetData(idst+i,j,k) = *simIso->lbmGetData(i,j,k); - } - } - } - - } // simIso - } - } - } -#endif // LBM_INCLUDE_TESTSOLVERS==1 + #if LBM_INCLUDE_TESTSOLVERS==1 + mrIsoExchange(); + #endif // LBM_INCLUDE_TESTSOLVERS==1 return; } @@ -471,12 +434,12 @@ int LbmFsgrSolver::initParticles() { partt->setEnd ( this->mvGeoEnd + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) ); partt->setSimStart( ntlVec3Gfx(0.0) ); - partt->setSimEnd ( ntlVec3Gfx(this->mSizex, this->mSizey, getForZMaxBnd(mMaxRefine)) ); + partt->setSimEnd ( ntlVec3Gfx(mSizex, mSizey, getForZMaxBnd(mMaxRefine)) ); while( (numgetNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) { LbmFloat x,y,z,t; - x = 1.0+(( (LbmFloat)(this->mSizex-3.) ) * (rand()/(RAND_MAX+1.0)) ); - y = 1.0+(( (LbmFloat)(this->mSizey-3.) ) * (rand()/(RAND_MAX+1.0)) ); + x = 1.0+(( (LbmFloat)(mSizex-3.) ) * (rand()/(RAND_MAX+1.0)) ); + y = 1.0+(( (LbmFloat)(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); @@ -485,8 +448,8 @@ int LbmFsgrSolver::initParticles() { k = 0; 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), CFNoNbFluid ) + //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFFluid) ) + //&& ( RFLAG(mMaxRefine, i,j,k, workSet)& CFNoNbFluid ) //if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) { if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid|CFUnused) ) { bool cellOk = true; @@ -531,8 +494,8 @@ int LbmFsgrSolver::initParticles() { #if LBM_INCLUDE_TESTSOLVERS==1 if(mUseTestdata) { const bool partDebug=false; - if(mpTest->mDebugvalue2!=0.0){ errMsg("LbmTestdata"," part init "<mDebugvalue2); } - if(mpTest->mDebugvalue2==-12.0){ + if(mpTest->mPartTestcase==0){ errMsg("LbmTestdata"," part init "<mPartTestcase); } + if(mpTest->mPartTestcase==-12){ const int lev = mMaxRefine; for(int i=5;i<15;i++) { LbmFloat x,y,z; @@ -546,7 +509,7 @@ int LbmFsgrSolver::initParticles() { if(partDebug) errMsg("PARTTT","SET "<getLast()->getPos() <<" s"<getLast()->getSize() ); } } - if(mpTest->mDebugvalue2==-11.0){ + if(mpTest->mPartTestcase==-11){ const int lev = mMaxRefine; for(int i=5;i<15;i++) { LbmFloat x,y,z; @@ -561,7 +524,7 @@ int LbmFsgrSolver::initParticles() { } } // place floats on rectangular region FLOAT_JITTER_BND - if(mpTest->mDebugvalue2==-10.0){ + if(mpTest->mPartTestcase==-10){ const int lev = mMaxRefine; const int sx = mLevel[lev].lSizex; const int sy = mLevel[lev].lSizey; @@ -594,6 +557,22 @@ int LbmFsgrSolver::initParticles() { return 0; } +// helper function for particle debugging +/*static string getParticleStatusString(int state) { + std::ostringstream out; + if(state&PART_DROP) out << "dropp "; + if(state&PART_TRACER) out << "tracr "; + if(state&PART_BUBBLE) out << "bubbl "; + if(state&PART_FLOAT) out << "float "; + if(state&PART_INTER) out << "inter "; + + if(state&PART_IN) out << "inn "; + if(state&PART_OUT) out << "out "; + if(state&PART_INACTIVE) out << "INACT "; + if(state&PART_OUTFLUID) out << "outfluid "; + return out.str(); +} // */ + #define P_CHANGETYPE(p, newtype) \ p->setLifeTime(0.); \ /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<getFlags())<<" to "<< (newtype) ); */ \ @@ -610,21 +589,22 @@ int LbmFsgrSolver::initParticles() { #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)) - -void LbmFsgrSolver::advanceParticles() { - int workSet = mLevel[mMaxRefine].setCurr; - LbmFloat vx=0.0,vy=0.0,vz=0.0; - //LbmFloat df[27]; //feq[27]; #define DEL_PART { \ /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<getType()<<" "); */ \ p->setActive( false ); \ continue; } +void LbmFsgrSolver::advanceParticles() { + const int level = mMaxRefine; + const int workSet = mLevel[level].setCurr; + LbmFloat vx=0.0,vy=0.0,vz=0.0; + //int debugOutCounter=0; // debug output counter + myTime_t parttstart = getTime(); const LbmFloat cellsize = this->mpParam->getCellSize(); const LbmFloat timestep = this->mpParam->getTimestep(); //const LbmFloat viscAir = 1.79 * 1e-5; // RW2L kin. viscosity, mu - const LbmFloat viscWater = 1.0 * 1e-6; // RW2L kin. viscosity, mu + //const LbmFloat viscWater = 1.0 * 1e-6; // RW2L kin. viscosity, mu const LbmFloat rhoAir = 1.2; // [kg m^-3] RW2L const LbmFloat rhoWater = 1000.0; // RW2L const LbmFloat minDropSize = 0.0005; // [m], = 2mm RW2L @@ -635,65 +615,56 @@ void LbmFsgrSolver::advanceParticles() { 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){ + const bool useff = (mFarFieldSize>1.2); // if(mpTest->mFarfMode>0){ - // TODO use timestep size + // TODO scale bubble/part damping dep. on timestep, also drop bnd rev damping const int cutval = mCutoff; // use full border!? - int actCnt=0; if(this->mStepCnt%50==49) { mpParticles->cleanup(); } for(vector::iterator pit= mpParticles->getParticlesBegin(); pit!= mpParticles->getParticlesEnd(); pit++) { - //errMsg("PIT"," pit"<<(*pit).getId()<<" pos:"<< (*pit).getPos()<<" status:"<getPos()<<" vel:"<< (*pit)->getVel()<<" status:"<getFlags()) <<" " <getStart()<<" "<getEnd() ); - //int flag = (*pit).getFlags(); + //if((*pit).getPos()[2]>10.) errMsg("PIT"," pit"<<(*pit).getId()<<" pos:"<< (*pit).getPos()<<" status:["<getLifeTime()<0.){ - if(p->getLifeTime()<-this->mSimulationTime) continue; - else p->setLifeTime(-mLevel[mMaxRefine].timestep); // zero for following update + if(p->getLifeTime() < -mSimulationTime) continue; + else p->setLifeTime(-mLevel[level].timestep); // zero for following update } int i,j,k; - p->setLifeTime(p->getLifeTime()+mLevel[mMaxRefine].timestep); + p->setLifeTime(p->getLifeTime()+mLevel[level].timestep); // nearest neighbor, particle positions don't include empty bounds ntlVec3Gfx pos = p->getPos(); - i= (int)(pos[0]+0.5); - j= (int)(pos[1]+0.5); - k= (int)(pos[2]+0.5); + i= (int)pos[0]; j= (int)pos[1]; k= (int)pos[2];// no offset necessary if(LBMDIM==2) { k = 0; } // only testdata handling, all for sws #if LBM_INCLUDE_TESTSOLVERS==1 - if(useff && (mpTest->mDebugvalue1>0.)) { + if(useff && (mpTest->mFarfMode>0)) { 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... + // in out tests if(p->getStatus()&PART_IN) { // IN - if( (ithis->mSizex-1-cutval)|| - (jthis->mSizey-1-cutval) - //||(kthis->mSizez-1-cutval) + if( (imSizex-1-cutval)|| + (jmSizey-1-cutval) + //||(kmSizez-1-cutval) ) { if(!useff) { DEL_PART; } else { p->setStatus(PART_OUT); - /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART; } } } else { // OUT rough check // check in again? - if( (i>=cutval)&&(i<=this->mSizex-1-cutval)&& - (j>=cutval)&&(j<=this->mSizey-1-cutval) - //&&(k>=cutval)&&(k<=this->mSizez-1-cutval) + if( (i>=cutval)&&(i<=mSizex-1-cutval)&& + (j>=cutval)&&(j<=mSizey-1-cutval) ) { p->setStatus(PART_IN); - /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART; } } - //p->setStatus(PART_OUT);// DEBUG always out! if( (p->getType()==PART_BUBBLE) || (p->getType()==PART_TRACER) ) { @@ -702,11 +673,11 @@ void LbmFsgrSolver::advanceParticles() { vx = vy = vz = 0.0; if(p->getStatus()&PART_IN) { // IN if(k>=cutval) { - if(k>this->mSizez-1-cutval) DEL_PART; + if(k>mSizez-1-cutval) DEL_PART; - if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFUnused) ) { + if( RFLAG(level, i,j,k, workSet)&(CFFluid|CFUnused) ) { // still ok - int partLev = mMaxRefine; + int partLev = level; int si=i, sj=j, sk=k; while(partLev>0 && RFLAG(partLev, si,sj,sk, workSet)&(CFUnused)) { partLev--; @@ -714,23 +685,22 @@ void LbmFsgrSolver::advanceParticles() { sj/=2; sk/=2; } - //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l); - LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet); - FORDF1{ - //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l); - LbmFloat cdf = RAC(ccel, l); - // TODO update below - //df[l] = cdf; - vx += (this->dfDvecX[l]*cdf); - vy += (this->dfDvecY[l]*cdf); - vz += (this->dfDvecZ[l]*cdf); - } - - // remove gravity influence - const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les - vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5; - vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5; - vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5; + // get velocity from fluid cell + if( RFLAG(partLev, si,sj,sk, workSet)&(CFFluid) ) { + LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet); + FORDF1{ + LbmFloat cdf = RAC(ccel, l); + // TODO update below + vx += (this->dfDvecX[l]*cdf); + vy += (this->dfDvecY[l]*cdf); + vz += (this->dfDvecZ[l]*cdf); + } + // remove gravity influence + const LbmFloat lesomega = mLevel[level].omega; // no les + vx -= mLevel[level].gravity[0] * lesomega*0.5; + vy -= mLevel[level].gravity[1] * lesomega*0.5; + vz -= mLevel[level].gravity[2] * lesomega*0.5; + } // fluid vel } else { // OUT // out of bounds, deactivate... @@ -753,39 +723,41 @@ void LbmFsgrSolver::advanceParticles() { ntlVec3Gfx v = p->getVel(); // dampen... if( (useff)&& (p->getType()==PART_BUBBLE) ) { // test rise - //O vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2]; - LbmFloat radius = p->getSize() * minDropSize; - LbmVec velPart = vec2L(p->getVel()) *cellsize/timestep; // L2RW, lattice velocity - LbmVec velWater = LbmVec(vx,vy,vz) *cellsize/timestep;// L2RW, fluid velocity - LbmVec velRel = velWater - velPart; - LbmFloat velRelNorm = norm(velRel); - // TODO calculate values in lattice units, compute CD?!??! - LbmFloat pvolume = rhoAir * 4.0/3.0 * M_PI* radius*radius*radius; // volume: 4/3 pi r^3 - //const LbmFloat cd = + if(mPartUsePhysModel) { + LbmFloat radius = p->getSize() * minDropSize; + LbmVec velPart = vec2L(p->getVel()) *cellsize/timestep; // L2RW, lattice velocity + LbmVec velWater = LbmVec(vx,vy,vz) *cellsize/timestep;// L2RW, fluid velocity + LbmVec velRel = velWater - velPart; + //LbmFloat velRelNorm = norm(velRel); + LbmFloat pvolume = rhoAir * 4.0/3.0 * M_PI* radius*radius*radius; // volume: 4/3 pi r^3 - LbmVec fb = -rwgrav* pvolume *rhoWater; - LbmVec fd = velRel*6.0*M_PI*radius* (1e-3); //viscWater; - LbmVec change = (fb+fd) *10.0*timestep *(timestep/cellsize); - //LbmVec change = (fb+fd) *timestep / (pvolume*rhoAir) *(timestep/cellsize); - //actCnt++; // should be after active test - if(actCnt<0) { - errMsg("PIT","BTEST1 vol="<getVel())) * 6.0*M_PI*radius* (1e-3); //viscWater; + LbmFloat w = 0.99; + vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[level].gravity[2]); + v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2); + p->setVel( v ); + } else { + // non phys, half old, half fluid, use slightly slower acc + v = v*0.5 + ntlVec3Gfx(vx,vy,vz)* 0.5-vec2G(mLevel[level].gravity)*0.5; + p->setVel( v * 0.99 ); } - - //v += change; - //v += ntlVec3Gfx(vx,vy,vz); - LbmVec fd2 = (LbmVec(vx,vy,vz)-vec2L(p->getVel())) * 6.0*M_PI*radius* (1e-3); //viscWater; - 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); + p->advanceVel(); + } else if(p->getType()==PART_TRACER) { v = ntlVec3Gfx(vx,vy,vz); - CellFlagType fflag = RFLAG(mMaxRefine, i,j,k, workSet); + CellFlagType fflag = RFLAG(level, i,j,k, workSet); if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true); } else { p->setInFluid(false); } @@ -801,16 +773,17 @@ void LbmFsgrSolver::advanceParticles() { } else { // move inwards along normal, make sure normal is valid first - const int lev = mMaxRefine; + // todo use class funcs! + const int lev = level; LbmFloat nx=0.,ny=0.,nz=0., nv1,nv2; bool nonorm = false; if(i<=0) { nx = -1.; nonorm = true; } - if(i>=this->mSizex-1) { nx = 1.; nonorm = true; } + if(i>=mSizex-1) { nx = 1.; nonorm = true; } if(j<=0) { ny = -1.; nonorm = true; } - if(j>=this->mSizey-1) { ny = 1.; nonorm = true; } + if(j>=mSizey-1) { ny = 1.; nonorm = true; } # if LBMDIM==3 if(k<=0) { nz = -1.; nonorm = true; } - if(k>=this->mSizez-1) { nz = 1.; nonorm = true; } + if(k>=mSizez-1) { nz = 1.; nonorm = true; } # endif // LBMDIM==3 if(!nonorm) { FFGET_NORM(nv1,dE); FFGET_NORM(nv2,dW); @@ -824,24 +797,21 @@ void LbmFsgrSolver::advanceParticles() { nz = 0.; # endif // LBMDIM==3 } else { - v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity); + v = p->getVel() + vec2G(mLevel[level].gravity); } - - p->advanceVec( (ntlVec3Gfx(nx,ny,nz)) * -0.1 ); // + vec2G(mLevel[mMaxRefine].gravity); - //if(normNoSqrt(v)advanceVec( (ntlVec3Gfx(nx,ny,nz)) * -0.1 ); // + vec2G(mLevel[level].gravity); } } p->setVel( v ); p->advanceVel(); - //errMsg("PPPP"," pos"<getPos()<<" "<getVel() ); } // drop handling else if(p->getType()==PART_DROP) { ntlVec3Gfx v = p->getVel(); // dampen... - if(1) { + if(mPartUsePhysModel) { LbmFloat radius = p->getSize() * minDropSize; LbmVec velPart = vec2L(p->getVel()) *cellsize /timestep; // * cellsize / timestep; // L2RW, lattice velocity LbmVec velRel = velAir - velPart; @@ -857,39 +827,98 @@ void LbmFsgrSolver::advanceParticles() { LbmVec fg = rwgrav * mb;// * (1.0-rhoAir/rhoWater); LbmVec fd = velRel* velRelNorm* cd*M_PI *rhoAir *0.5 *radius*radius; LbmVec change = (fg+ fd ) *timestep / mb *(timestep/cellsize); - - //actCnt++; // should be after active test - if(actCnt<0) { - errMsg("\nPIT","NTEST1 mb="<0) { errMsg("\nPIT","NTEST1 mb="<setVel( v*0.9 ); // restdamping + } else { + p->addToVel( vec2G(mLevel[level].gravity) ); + } } // OLD p->advanceVel(); if(p->getStatus()&PART_IN) { // IN if(kmSizez-1-cutval){ - //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter)) { - if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) { + if(k<=mSizez-1-cutval){ + CellFlagType pflag = RFLAG(level, i,j,k, workSet); + //errMsg("PIT move"," at "<getSize(); + const LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*size); + bool orgcellok = false; + if( RFLAG(level, oi,oj,ok, workSet) & (CFInter) ){ + orgcellok = true; + } else { + // search upward for interface + oi=i; oj=j; ok=k; + for(int kk=0; kk<5 && ok<=mSizez-1; kk++) { + ok++; + if( RFLAG(level, oi,oj,ok, workSet) & (CFInter) ){ + kk = 5; orgcellok = true; + } + } + } + + //errMsg("PTIMPULSE"," new v"< 0.166*0.166) { + v *= 1./sqrtf(vlensqr)*0.166; + } + // compute cell velocity + LbmFloat cellVelSqr = 0.; + LbmFloat *tcel = RACPNT(level, oi,oj,ok, workSet); + LbmFloat velUx=0., velUy=0., velUz=0.; + FORDF0 { + velUx += (this->dfDvecX[l]*RAC(tcel,l)); + velUy += (this->dfDvecY[l]*RAC(tcel,l)); + velUz += (this->dfDvecZ[l]*RAC(tcel,l)); + } + // add impulse + /*cellVelSqr = velUx*velUx+ velUy*velUy+ velUz*velUz; + //errMsg("PTIMPULSE"," new v"<getType()==PART_FLOAT) { - // test - delte on boundary!? - //if( (i<=cutval)||(i>=this->mSizex-1-cutval)|| (j<=cutval)||(j>=this->mSizey-1-cutval)) { DEL_PART; } // DEBUG TEST - -# if LBM_INCLUDE_TESTSOLVERS==1 - /*LbmFloat prob = 1.0; - // vanishing - prob = (rand()/(RAND_MAX+1.0)); - // 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(probthis->mSizez*3/4) { if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;} - // */ - -# endif // LBM_INCLUDE_TESTSOLVERS==1 if(p->getStatus()&PART_IN) { // IN - //if((kthis->mSizez-1-cutval)) DEL_PART; if(kthis->mSizez-1-cutval) { P_CHANGETYPE(p, PART_DROP ); continue; } // create new drop - //ntlVec3Gfx v = getVelocityAt(i,j,k); + // not valid for mass... vx = vy = vz = 0.0; // define from particletracer.h @@ -973,11 +974,10 @@ void LbmFsgrSolver::advanceParticles() { int ccnt=0; for(int kk=0;kkdfDvecX[l]*cdf); vy += (this->dfDvecY[l]*cdf); vz += (this->dfDvecZ[l]*cdf); @@ -995,55 +995,26 @@ void LbmFsgrSolver::advanceParticles() { vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5); //bool delfloat = false; - if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) { + if( ( RFLAG(level, i,j,k, workSet)& (CFFluid|CFUnused) ) ) { // in fluid cell - /*if((1) && (kmSizez-3) && - ( - 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; - // keep below obstacles - if((delfloat) && (TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFBnd )) ) { - delfloat=false; vz=0.0; - } // */ - vz = p->getVel()[2]-1.0*mLevel[mMaxRefine].gravity[2]; // simply rise... + vz = p->getVel()[2]-1.0*mLevel[level].gravity[2]; // simply rise... if(vz<0.) vz=0.; - } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFBnd ) ) { - //vz = p->getVel()[2]-0.2; // force downwards movement, move below obstacle... - vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall... - if(vz>0.) vz=0.; - } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) { + } else if( ( RFLAG(level, i,j,k, workSet)& CFBnd ) ) { + // force downwards movement, move below obstacle... + //vz = p->getVel()[2]+1.0*mLevel[level].gravity[2]; // fall... + //if(vz>0.) vz=0.; + DEL_PART; + } else if( ( RFLAG(level, 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 - vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall... - if(vz>0.) vz=0.; - // 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; // */ - //P_CHANGETYPE(p, PART_DROP ); continue; // create new drop + } else { // all else... + if( ( RFLAG(level, i,j,k-1, workSet)& (CFFluid|CFInter) ) ) { + vz = p->getVel()[2]+2.0*mLevel[level].gravity[2]; // fall... + if(vz>0.) vz=0.; } + else { DEL_PART; } } - //if(delfloat) DEL_PART; - /* - // move down from empty - else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { - vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2]; - if(vz>0.0) vz=0.0; - //DEL_PART; // ???? - } else { DEL_PART; } // */ - //vz = 0.0; // DEBUG 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:"<handleParticle(p, i,j,k); } @@ -1051,19 +1022,17 @@ void LbmFsgrSolver::advanceParticles() { #else // LBM_INCLUDE_TESTSOLVERS==1 DEL_PART; #endif // LBM_INCLUDE_TESTSOLVERS==1 - //errMsg("PIT","OUT pit "<< (*pit).getPos()<<" status:"<getLifeTime()<3.*mLevel[mMaxRefine].timestep)) { + if((0) && (useff) && (p->getLifeTime()<3.*mLevel[level].timestep)) { // use half butoff border 1/8 - int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5); + int maxdw = (int)(mLevel[level].lSizex*0.125*0.5); if(maxdw<3) maxdw=3; - if((j>=0)&&(j<=this->mSizey-1)) { + if((j>=0)&&(j<=mSizey-1)) { if(ABS(i-( cutval))advance( FLOAT_JITTBNDRAND( ABS(i-( cutval))), 0.,0.); } - if(ABS(i-(this->mSizex-1-cutval))advance( FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); } + if(ABS(i-(mSizex-1-cutval))advance( FLOAT_JITTBNDRAND( ABS(i-(mSizex-1-cutval))), 0.,0.); } } - //if( (ithis->mSizex-1-cutval)|| //(jthis->mSizey-1-cutval) } } // PART_FLOAT @@ -1080,9 +1049,8 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt int workSet = mLevel[mMaxRefine].setCurr; std::ostringstream name; - // debug - raw dump of ffrac values - if(getenv("ELBEEM_RAWDEBUGDUMP")) { - //name <<"fill_" << this->mStepCnt <<".dump"; + // debug - raw dump of ffrac values, as text! + if(mDumpRawText) { name << outfilename<< frameNrStr <<".dump"; FILE *file = fopen(name.str().c_str(),"w"); if(file) { @@ -1108,31 +1076,156 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt } // file } // */ - if(getenv("ELBEEM_BINDEBUGDUMP")) { - name << outfilename<< frameNrStr <<".bdump"; - FILE *file = fopen(name.str().c_str(),"w"); - if(file) { - for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { - for(int j=0;j1.) val=1.; + + if(mDumpRawBinary) { + if(!mDumpRawBinaryZip) { + // unzipped, only fill + name << outfilename<< frameNrStr <<".bdump"; + FILE *file = fopen(name.str().c_str(),"w"); + if(file) { + for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { + for(int j=0;j1.) val=1.; + } + if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; + fwrite( &val, sizeof(val), 1, file); // binary } - if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; - fwrite( &val, sizeof(val), 1, file); // binary } } - } - fclose(file); - } // file - } + fclose(file); + } // file + } // unzipped + else { + // zipped, use iso values + prepareVisualization(); + name << outfilename<< frameNrStr <<".bdump.gz"; + gzFile gzf = gzopen(name.str().c_str(),"wb9"); + if(gzf) { + // write size + int s; + s=mSizex; gzwrite(gzf, &s, sizeof(s)); + s=mSizey; gzwrite(gzf, &s, sizeof(s)); + s=mSizez; gzwrite(gzf, &s, sizeof(s)); + + // write isovalues + for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { + for(int j=0;jlbmGetData( i,j,k ); + gzwrite(gzf, &val, sizeof(val)); + } + } + } + gzclose(gzf); + } // gzf + } // zip + } // bin dump dumptype = 0; frameNr = 0; // get rid of warning } +/*! move a particle at a boundary */ +void LbmFsgrSolver::handleObstacleParticle(ParticleObject *p) { + //if(normNoSqrt(v)<=0.) continue; // skip stuck + /* + p->setVel( v * -1. ); // revert + p->advanceVel(); // move back twice... + if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFBndNoslip)) { + p->setVel( v * -0.5 ); // revert & dampen + } + p->advanceVel(); + // */ + // TODO mark/remove stuck parts!? + + const int level = mMaxRefine; + const int workSet = mLevel[level].setCurr; + LbmVec v = vec2L( p->getVel() ); + if(normNoSqrt(v)<=0.) { + p->setVel(vec2G(mLevel[level].gravity)); + } + + CellFlagType pflag = CFBnd; + ntlVec3Gfx posOrg(p->getPos()); + ntlVec3Gfx npos(0.); + int ni=1,nj=1,nk=1; + int tries = 0; + + // try to undo movement + p->advanceVec( (p->getVel()-vec2G(mLevel[level].gravity)) * -2.); + + npos = p->getPos(); ni= (int)npos[0]; + nj= (int)npos[1]; nk= (int)npos[2]; + if(LBMDIM==2) { nk = 0; } + //errMsg("BOUNDCPAR"," t"<setActive( false ); + return; + } + pflag = RFLAG(level, ni,nj,nk, workSet); + + // try to force particle out of boundary + bool haveNorm = false; + LbmVec bnormal; + if(pflag&CFBnd) { + npos = posOrg; ni= (int)npos[0]; + nj= (int)npos[1]; nk= (int)npos[2]; + if(LBMDIM==2) { nk = 0; } + + computeObstacleSurfaceNormalAcc(ni,nj,nk, &bnormal[0]); + haveNorm = true; + normalize(bnormal); + bnormal *= 0.25; + + tries = 1; + while(pflag&CFBnd && tries<=5) { + // use increasing step sizes + p->advanceVec( vec2G( bnormal *0.5 *(gfxReal)tries ) ); + npos = p->getPos(); + ni= (int)npos[0]; + nj= (int)npos[1]; + nk= (int)npos[2]; + + // delete out of domain + if(!checkDomainBounds(level,ni,nj,nk)) { + //errMsg("BOUNDCPAR"," DEL! "); + p->setActive( false ); + return; + } + pflag = RFLAG(level, ni,nj,nk, workSet); + tries++; + } + + // really stuck, delete... + if(pflag&CFBnd) { + p->setActive( false ); + return; + } + } + + // not in bound anymore! + if(!haveNorm) { + CellFlagType *bflag = &RFLAG(level, ni,nj,nk, workSet); + LbmFloat *bcell = RACPNT(level, ni,nj,nk, workSet); + computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]); + } + normalize(bnormal); + LbmVec normComp = bnormal * dot(vec2L(v),bnormal); + //errMsg("BOUNDCPAR","bnormal"<setVel(vec2G(v)); + p->advanceVel(); +} + /*****************************************************************************/ /*! internal quick print function (for debugging) */ /*****************************************************************************/ @@ -1320,6 +1413,13 @@ ntlVec3Gfx LbmFsgrSolver::getCellSize ( CellIdentifierInterface* basecid) LbmFloat LbmFsgrSolver::getCellDensity ( CellIdentifierInterface* basecid,int set) { stdCellId *cid = convertBaseCidToStdCid(basecid); + // skip non-fluid cells + if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&(CFFluid|CFInter)) { + // ok go on... + } else { + return 0.; + } + LbmFloat rho = 0.0; FORDF0 { rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l); } // ORG return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep) +1.0; // ORG @@ -1744,7 +1844,7 @@ void LbmFsgrSolver::lbmDebugDisplay(int dispset) { glDisable( GL_LIGHTING ); // dont light lines #if LBM_INCLUDE_TESTSOLVERS==1 - if((!mUseTestdata)|| (mUseTestdata)&&(mpTest->mDebugvalue1<=0.0)) { + if((!mUseTestdata)|| (mUseTestdata)&&(mpTest->mFarfMode<=0)) { #endif // LBM_INCLUDE_TESTSOLVERS==1 LbmFsgrSolver::CellIdentifier cid = this->getFirstCell(); diff --git a/intern/elbeem/intern/utilities.cpp b/intern/elbeem/intern/utilities.cpp index 1aee43fe8bd..ad827efa2fc 100644 --- a/intern/elbeem/intern/utilities.cpp +++ b/intern/elbeem/intern/utilities.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Global C style utility funcions * @@ -28,6 +28,7 @@ #include #endif #endif // NOPNG +#include // global debug level #ifdef DEBUG @@ -64,16 +65,28 @@ char* getElbeemErrorString(void) { return gElbeemErrorString; } myTime_t globalIntervalTime = 0; //! color output setting for messages (0==off, else on) #ifdef WIN32 -int globalColorSetting = 0; +// switch off first call +#define DEF_globalColorSetting -1 #else // WIN32 -int globalColorSetting = 1; +// linux etc., on by default +#define DEF_globalColorSetting 1 #endif // WIN32 +int globalColorSetting = DEF_globalColorSetting; // linux etc., on by default int globalFirstEnvCheck = 0; +void resetGlobalColorSetting() { globalColorSetting = DEF_globalColorSetting; } // global string for formatting vector output, TODO test!? char *globVecFormatStr = "V[%f,%f,%f]"; +// global mp on/off switch +bool glob_mpactive = false; +// global access to mpi index, for debugging (e.g. in utilities.cpp) +int glob_mpnum = -1; +int glob_mpindex = -1; +int glob_mppn = -1; + + //----------------------------------------------------------------------------- // helper function that converts a string to integer, // and returns an alternative value if the conversion fails @@ -126,7 +139,7 @@ int writePng(const char *fileName, unsigned char **rowsp, int w, int h) //FILE *fp = fopen(fileName, "wb"); FILE *fp = NULL; - char *doing = "open for writing"; + string doing = "open for writing"; if (!(fp = fopen(fileName, "wb"))) goto fail; if(!png_ptr) { @@ -162,6 +175,40 @@ fail: if(png_ptr || info_ptr) png_destroy_write_struct(&png_ptr, &info_ptr); return -1; } +#else // NOPNG +// fallback - write ppm +int writePng(const char *fileName, unsigned char **rowsp, int w, int h) +{ + gzFile gzf; + string filentemp(fileName); + // remove suffix + if((filentemp.length()>4) && (filentemp[filentemp.length()-4]=='.')) { + filentemp[filentemp.length()-4] = '\0'; + } + std::ostringstream filennew; + filennew << filentemp.c_str(); + filennew << ".ppm.gz"; + + gzf = gzopen(filennew.str().c_str(), "wb9"); + if(!gzf) goto fail; + + gzprintf(gzf,"P6\n%d %d\n255\n",w,h); + // output binary pixels + for(int j=0;j0) { myTime_t currTime = getTime(); if((currTime - globalIntervalTime)>interval) { @@ -254,14 +315,16 @@ void messageOutputFunc(string from, int id, string msg, myTime_t interval) { } // colors off? - if((globalColorSetting==0) || (id==DM_FATAL) ){ + if( (globalColorSetting == -1) || // off for e.g. win32 + ((globalColorSetting==1) && ((id==DM_FATAL)||( getenv("ELBEEM_NOCOLOROUT") )) ) + ) { // only reset once col_std = col_black = col_dark_gray = col_bright_gray = col_red = col_bright_red = col_green = col_bright_green = col_bright_yellow = col_yellow = col_cyan = col_bright_cyan = col_purple = col_bright_purple = col_neutral = ""; - globalColorSetting=1; + globalColorSetting = 0; } std::ostringstream sout; @@ -288,6 +351,9 @@ void messageOutputFunc(string from, int id, string msg, myTime_t interval) { case DM_FATAL: sout << col_red << " fatal("<0 anyway, so write to file... - // TODO generate some reasonable path? - FILE *logf = fopen("elbeem_debug_log.txt","a+"); - // dont complain anymore here... - if(logf) { - fprintf(logf, "%s",sout.str().c_str() ); - fclose(logf); + // determine output - file==1/stdout==0 / globstr==2 + char filen[256]; + strcpy(filen,"debug_unini.txt"); + int fileout = false; +#if ELBEEM_MPI==1 + std::ostringstream mpin; + if(glob_mpindex>=0) { + mpin << "elbeem_log_"<< glob_mpindex <<".txt"; + } else { + mpin << "elbeem_log_ini.txt"; } -#else // WIN32 - fprintf(stdout, "%s",sout.str().c_str() ); - if(id!=DM_DIRECT) fflush(stdout); + fileout = 1; + strncpy(filen, mpin.str().c_str(),255); filen[255]='\0'; +#else + strncpy(filen, "elbeem_debug_log.txt",255); +#endif + +#ifdef WIN32 + // windows causes trouble with direct output + fileout = 1; #endif // WIN32 + +#if PARALLEL==1 + fileout = 2;// buffer out, switch off again... + if(globOutstrForce) fileout=1; +#endif + if(getenv("ELBEEM_FORCESTDOUT")) { + fileout = 0;// always direct out + } + //fprintf(stdout,"out deb %d, %d, '%s',l%d \n",globOutstrForce,fileout, filen, globOutstr.str().size() ); + +#if PARALLEL==1 +#pragma omp critical +#endif // PARALLEL==1 + { + if(fileout==1) { + // debug level is >0 anyway, so write to file... + FILE *logf = fopen(filen,"a+"); + // dont complain anymore here... + if(logf) { + if(globOutstrForce) { + fprintf(logf, "%s",globOutstr.str().c_str() ); + globOutstr.str(""); // reset + } + fprintf(logf, "%s",sout.str().c_str() ); + fclose(logf); + } + } else if(fileout==2) { + globOutstr << sout.str(); + } else { + // normal stdout output + fprintf(stdout, "%s",sout.str().c_str() ); + if(id!=DM_DIRECT) fflush(stdout); + } + } // omp crit } // helper functions from external program using elbeem lib (e.g. Blender) @@ -323,14 +431,19 @@ void messageOutputFunc(string from, int id, string msg, myTime_t interval) { extern "C" void elbeemCheckDebugEnv(void) { const char *strEnvName = "BLENDER_ELBEEMDEBUG"; + const char *strEnvName2 = "ELBEEM_DEBUGLEVEL"; if(globalFirstEnvCheck) return; if(getenv(strEnvName)) { gDebugLevel = atoi(getenv(strEnvName)); - if(gDebugLevel< 0) gDebugLevel = 0; - if(gDebugLevel>10) gDebugLevel = 0; // only use valid values if(gDebugLevel>0) debMsgStd("performElbeemSimulation",DM_NOTIFY,"Using envvar '"<0) debMsgStd("performElbeemSimulation",DM_NOTIFY,"Using envvar '"<10) gDebugLevel = 0; // only use valid values globalFirstEnvCheck = 1; } @@ -367,7 +480,8 @@ double elbeemEstimateMemreq(int res, double memreq = -1.0; string memreqStr(""); - calculateMemreqEstimate(resx,resy,resz, refine, &memreq, &memreqStr ); + // ignore farfield for now... + calculateMemreqEstimate(resx,resy,resz, refine, 0., &memreq, &memreqStr ); if(retstr) { // copy at max. 32 characters diff --git a/intern/elbeem/intern/utilities.h b/intern/elbeem/intern/utilities.h index 57a956c64d9..4d887c3f99b 100644 --- a/intern/elbeem/intern/utilities.h +++ b/intern/elbeem/intern/utilities.h @@ -1,7 +1,7 @@ /****************************************************************************** * * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method - * Copyright 2003,2004 Nils Thuerey + * Copyright 2003-2006 Nils Thuerey * * Global C style utility funcions * @@ -115,6 +115,8 @@ string getTimeString(myTime_t usecs); //! helper to check if a bounding box was specified in the right way bool checkBoundingBox(ntlVec3Gfx s, ntlVec3Gfx e, string checker); +//! reset color output for elbeem init +void resetGlobalColorSetting(); /*! print some vector from 3 values e.g. for ux,uy,uz */ @@ -144,10 +146,8 @@ bool checkBoundingBox(ntlVec3Gfx s, ntlVec3Gfx e, string checker); PRINT_VEC( (mpV[(t).getPoints()[2]][0]),(mpV[(t).getPoints()[2]][1]),(mpV[(t).getPoints()[2]][2]) )<<" } " -#ifndef NOPNG // write png image int writePng(const char *fileName, unsigned char **rowsp, int w, int h); -#endif // NOPNG /* some useful templated functions * may require some operators for the classes diff --git a/source/blender/blenkernel/intern/DerivedMesh.c b/source/blender/blenkernel/intern/DerivedMesh.c index 3ccdbe130a6..4dacfdc85d8 100644 --- a/source/blender/blenkernel/intern/DerivedMesh.c +++ b/source/blender/blenkernel/intern/DerivedMesh.c @@ -917,6 +917,7 @@ static DerivedMesh *getMeshDerivedMesh(Mesh *me, Object *ob, float (*vertCos)[3] if((ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && (ob->fluidsimSettings->type & OB_FLUIDSIM_DOMAIN)&& (ob->fluidsimSettings->meshSurface) && + (1) && (!give_parteff(ob)) && // doesnt work together with particle systems! (me->totvert == ((Mesh *)(ob->fluidsimSettings->meshSurface))->totvert) ) { // dont recompute for fluidsim mesh, use from readBobjgz // TODO? check for modifiers!? @@ -2816,7 +2817,8 @@ static void mesh_calc_modifiers(Object *ob, float (*inputVertexCos)[3], * domain objects */ if((G.obedit!=ob) && !needMapping) { - if (ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) { + if((ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && + (1) && (!give_parteff(ob)) ) { // doesnt work together with particle systems! if(ob->fluidsimSettings->type & OB_FLUIDSIM_DOMAIN) { loadFluidsimMesh(ob,useRenderParams); fluidsimMeshUsed = 1; diff --git a/source/blender/blenkernel/intern/effect.c b/source/blender/blenkernel/intern/effect.c index a22bdad628e..08ea3b9abe1 100644 --- a/source/blender/blenkernel/intern/effect.c +++ b/source/blender/blenkernel/intern/effect.c @@ -1685,11 +1685,11 @@ void build_particle_system(Object *ob) if(paf->keys) MEM_freeN(paf->keys); /* free as early as possible, for returns */ paf->keys= NULL; - printf("build particles\n"); + //printf("build particles\n"); - /* fluid sim particle import handling, actual loading */ + /* fluid sim particle import handling, actual loading of particles from file */ #ifndef DISABLE_ELBEEM - if( (1) && (ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && + if( (1) && (ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) && // broken, disabled for now! (ob->fluidsimSettings) && (ob->fluidsimSettings->type == OB_FLUIDSIM_PARTICLE)) { char *suffix = "fluidsurface_particles_#"; @@ -1702,7 +1702,7 @@ void build_particle_system(Object *ob) float vel[3]; if(ob==G.obedit) { // off... - paf->totpart = 1; + paf->totpart = 0; // 1 or 0? return; } @@ -1715,8 +1715,10 @@ void build_particle_system(Object *ob) gzf = gzopen(filename, "rb"); if (!gzf) { snprintf(debugStrBuffer,256,"readFsPartData::error - Unable to open file for reading '%s' \n", filename); +fprintf(stderr,"readFsPartData::error - Unable to open file for reading '%s' \n", filename); // DEBUG! //elbeemDebugOut(debugStrBuffer); - paf->totpart = 1; + //fprintf(stderr,"NO PARTOP!\n"); + paf->totpart = 0; return; } @@ -1753,7 +1755,7 @@ void build_particle_system(Object *ob) // convert range of 1.0-10.0 to shorts 1000-10000) shsize = (short)(convertSize*1000.0); pa->rt = shsize; - //if(a<200) fprintf(stderr,"SREAD %f %d %d \n",convertSize,shsize,pa->rt); + //if(a<200) fprintf(stderr,"SREAD, %d/%d: %f %d %d \n",a,totpart, convertSize,shsize,pa->rt); for(j=0; j<3; j++) { float wrf; diff --git a/source/blender/makesdna/DNA_object_fluidsim.h b/source/blender/makesdna/DNA_object_fluidsim.h index 5ddfc22045d..91c25d2ec1d 100644 --- a/source/blender/makesdna/DNA_object_fluidsim.h +++ b/source/blender/makesdna/DNA_object_fluidsim.h @@ -109,6 +109,9 @@ typedef struct FluidsimSettings { float generateParticles; /* smooth fluid surface? */ float surfaceSmoothing; + /* number of surface subdivisions*/ + int surfaceSubdivs; + int unusedDNADummy; /* particle display - size scaling, and alpha influence */ float particleInfSize, particleInfAlpha; diff --git a/source/blender/src/buttons_object.c b/source/blender/src/buttons_object.c index 65cf065bd66..8f05e20af0c 100644 --- a/source/blender/src/buttons_object.c +++ b/source/blender/src/buttons_object.c @@ -1608,21 +1608,24 @@ void do_object_panels(unsigned short event) break; case B_FLUIDSIM_MAKEPART: ob= OBACT; - { + if(1) { PartEff *paf= NULL; /* prepare fluidsim particle display */ // simplified delete effect, create new - recalc some particles... if(ob==NULL || ob->type!=OB_MESH) break; ob->fluidsimSettings->type = 0; // reset type, and init particle system once normally - eff= ob->effect.first; - //if((eff) && (eff->flag & SELECT)) { BLI_remlink(&ob->effect, eff); free_effect(eff); } - if(!eff){ copy_act_effect(ob); DAG_scene_sort(G.scene); } + paf= give_parteff(ob); + if( (BLI_countlist(&ob->effect)totpart = 1000; paf->sta = paf->end = 1.0; // generate some particles... - build_particle_system(ob); - - ob->fluidsimSettings->type = OB_FLUIDSIM_PARTICLE; + if(paf) { + paf->totpart = 1000; paf->sta = paf->end = 1.0; // generate some particles... + build_particle_system(ob); + ob->fluidsimSettings->type = OB_FLUIDSIM_PARTICLE; + } } allqueue(REDRAWVIEW3D, 0); allqueue(REDRAWBUTSOBJECT, 0); @@ -1977,6 +1980,14 @@ void do_effects_panels(unsigned short event) break; case B_NEWEFFECT: if(ob) { + if(ob->fluidsimFlag & OB_FLUIDSIM_ENABLE) { + // NT particles and fluid meshes currently dont work together -> switch off beforehand + if(ob->fluidsimSettings->type == OB_FLUIDSIM_DOMAIN) { + pupmenu("Fluidsim Particle Error%t|Please disable the fluidsim domain before activating particles.%x0"); + break; + //ob->fluidsimFlag = 0; ob->fluidsimSettings->type = 0; allqueue(REDRAWVIEW3D, 0); + } + } if (BLI_countlist(&ob->effect)==MAX_EFFECT) error("Unable to add: effect limit reached"); else @@ -2610,8 +2621,8 @@ static void object_panel_fluidsim(Object *ob) uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Obstacle", 230, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_OBSTACLE,20.0, 3.0, "Object is a fixed obstacle."); yline -= lineHeight; - uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Inflow", 90, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_INFLOW, 20.0, 4.0, "Object adds fluid to the simulation."); - uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Outflow", 160, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_OUTFLOW, 20.0, 5.0, "Object removes fluid from the simulation."); + uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Inflow", 90, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_INFLOW, 20.0, 4.0, "Object adds fluid to the simulation."); + uiDefButS(block, ROW, REDRAWBUTSOBJECT ,"Outflow", 160, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_OUTFLOW, 20.0, 5.0, "Object removes fluid from the simulation."); uiDefButS(block, ROW, B_FLUIDSIM_MAKEPART ,"Particle", 230, yline, 70,objHeight, &fss->type, 15.0, OB_FLUIDSIM_PARTICLE,20.0, 3.0, "Object is made a particle system to display particles generated by a fluidsim domain object."); uiBlockEndAlign(block); yline -= lineHeight; @@ -2689,7 +2700,7 @@ static void object_panel_fluidsim(Object *ob) uiDefButS(block, MENU, REDRAWVIEW3D, "Viscosity%t|Manual %x1|Water %x2|Oil %x3|Honey %x4", 0,yline, 90,objHeight, &fss->viscosityMode, 0, 0, 0, 0, "Set viscosity of the fluid to a preset value, or use manual input."); if(fss->viscosityMode==1) { - uiDefButF(block, NUM, B_DIFF, "Value:", 90, yline, 105,objHeight, &fss->viscosityValue, 0.0, 1.0, 10, 0, "Viscosity setting, value that is multiplied by 10 to the power of (exponent*-1)."); + uiDefButF(block, NUM, B_DIFF, "Value:", 90, yline, 105,objHeight, &fss->viscosityValue, 0.0, 10.0, 10, 0, "Viscosity setting, value that is multiplied by 10 to the power of (exponent*-1)."); uiDefButS(block, NUM, B_DIFF, "Neg-Exp.:", 195, yline, 105,objHeight, &fss->viscosityExponent, 0, 10, 10, 0, "Negative exponent for the viscosity value (to simplify entering small values e.g. 5*10^-6."); uiBlockEndAlign(block); } else { @@ -2715,8 +2726,8 @@ static void object_panel_fluidsim(Object *ob) } else if(fss->show_advancedoptions == 2) { // copied from obstacle... //yline -= lineHeight + 5; - uiDefBut(block, LABEL, 0, "Domain boundary type settings:", 0,yline,300,objHeight, NULL, 0.0, 0, 0, 0, ""); - yline -= lineHeight; + //uiDefBut(block, LABEL, 0, "Domain boundary type settings:", 0,yline,300,objHeight, NULL, 0.0, 0, 0, 0, ""); + //yline -= lineHeight; uiBlockBeginAlign(block); // domain 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."); @@ -2725,11 +2736,13 @@ static void object_panel_fluidsim(Object *ob) uiBlockEndAlign(block); yline -= lineHeight; - uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); if(fss->typeFlags&OB_FSBND_PARTSLIP) { + uiDefBut(block, LABEL, 0, "PartSlipValue:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); uiDefButF(block, NUM, B_DIFF, "", 200, yline,100,objHeight, &fss->partSlipValue, 0.0, 1.0, 10,0, "."); - } else { uiDefBut(block, LABEL, 0, "-", 200,yline,100,objHeight, NULL, 0.0, 0, 0, 0, ""); } - yline -= lineHeight; + yline -= lineHeight; + } else { + //uiDefBut(block, LABEL, 0, "-", 200,yline,100,objHeight, NULL, 0.0, 0, 0, 0, ""); + } // copied from obstacle... uiDefBut(block, LABEL, 0, "Tracer Particles:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); @@ -2738,6 +2751,9 @@ static void object_panel_fluidsim(Object *ob) 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, "Surface Subdiv:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); + uiDefButI(block, NUM, B_DIFF, "", 200, yline,100,objHeight, &fss->surfaceSubdivs, 0.0, 5.0, 10,0, "Number of isosurface subdivisions. This is necessary for the inclusion of particles into the surface generation. Warning - can lead to longer computation times!"); + yline -= lineHeight; 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)."); @@ -2777,6 +2793,12 @@ static void object_panel_fluidsim(Object *ob) yline -= lineHeight; } else { } + + // domainNovecgen "misused" here + uiDefBut(block, LABEL, 0, "Animated Mesh:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); + uiDefButBitC(block, TOG, 1, REDRAWBUTSOBJECT, "Export", 200, yline,100,objHeight, &fss->domainNovecgen, 0, 0, 0, 0, "Export this mesh as an animated one. Slower, only use if really necessary (e.g. armatures or parented objects), animated pos/rot/scale IPOs do not require it."); + yline -= lineHeight; + } // fluid inflow else if( (fss->type == OB_FLUIDSIM_OUTFLOW) ) { yline -= lineHeight + 5; @@ -2788,6 +2810,11 @@ static void object_panel_fluidsim(Object *ob) uiBlockEndAlign(block); yline -= lineHeight; + // domainNovecgen "misused" here + uiDefBut(block, LABEL, 0, "Animated Mesh:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); + uiDefButBitC(block, TOG, 1, REDRAWBUTSOBJECT, "Export", 200, yline,100,objHeight, &fss->domainNovecgen, 0, 0, 0, 0, "Export this mesh as an animated one. Slower, only use if really necessary (e.g. armatures or parented objects), animated pos/rot/scale IPOs do not require it."); + yline -= lineHeight; + //uiDefBut(block, LABEL, 0, "No additional settings as of now...", 0,yline,300,objHeight, NULL, 0.0, 0, 0, 0, ""); } else if( (fss->type == OB_FLUIDSIM_OBSTACLE) ) { @@ -2801,7 +2828,7 @@ static void object_panel_fluidsim(Object *ob) yline -= lineHeight; uiBlockBeginAlign(block); // obstacle - 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 ,"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); @@ -2809,7 +2836,7 @@ static void object_panel_fluidsim(Object *ob) // domainNovecgen "misused" here uiDefBut(block, LABEL, 0, "Animated Mesh:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); - uiDefButBitC(block, TOG, 1, REDRAWBUTSOBJECT, "Export", 200, yline,100,objHeight, &fss->domainNovecgen, 0, 0, 0, 0, "Export this mesh as an animated one. Slower, only use if really necessary (e.g. armatures), animated pos/rot/scale IPOs do not require it."); + uiDefButBitC(block, TOG, 1, REDRAWBUTSOBJECT, "Export", 200, yline,100,objHeight, &fss->domainNovecgen, 0, 0, 0, 0, "Export this mesh as an animated one. Slower, only use if really necessary (e.g. armatures or parented objects), animated pos/rot/scale IPOs do not require it."); yline -= lineHeight; uiDefBut(block, LABEL, 0, "PartSlip Amount:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); @@ -2818,12 +2845,25 @@ static void object_panel_fluidsim(Object *ob) } else { uiDefBut(block, LABEL, 0, "-", 200,yline,100,objHeight, NULL, 0.0, 0, 0, 0, ""); } yline -= lineHeight; + // generateParticles "misused" here + uiDefBut(block, LABEL, 0, "Impact Factor:", 0,yline,200,objHeight, NULL, 0.0, 0, 0, 0, ""); + uiDefButF(block, NUM, B_DIFF, "", 200, yline,100,objHeight, &fss->surfaceSmoothing, -2.0, 10.0, 10,0, "This is an unphsyical value for moving objects - it control the impact an obstacle has on the fluid, =0 behaves a bit like outflow (deleting fluid), =1 is default, while >1 results in high forces. Can be used to tweak total mass."); + yline -= lineHeight; + yline -= lineHeight; // obstacle } else if(fss->type == OB_FLUIDSIM_PARTICLE) { + + //fss->type == 0; // off, broken... + if(1) { + // limited selection, old fixed: fss->typeFlags = (1<<5)|(1<<1); +# define PARTBUT_WIDTH (300/3) + uiDefButBitS(block, TOG, (1<<2) , REDRAWBUTSOBJECT, "Drops", 0*PARTBUT_WIDTH, yline, PARTBUT_WIDTH,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Show drop particles."); + uiDefButBitS(block, TOG, (1<<4) , REDRAWBUTSOBJECT, "Floats", 1*PARTBUT_WIDTH, yline, PARTBUT_WIDTH,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Show floating foam particles."); + uiDefButBitS(block, TOG, (1<<5) , REDRAWBUTSOBJECT, "Tracer", 2*PARTBUT_WIDTH, yline, PARTBUT_WIDTH,objHeight, &fss->typeFlags, 0, 0, 0, 0, "Show tracer particles."); + yline -= lineHeight; +# undef PARTBUT_WIDTH - // fixed for now - fss->typeFlags = (1<<5)|(1<<1); uiDefBut(block, LABEL, 0, "Size Influence:", 0,yline,150,objHeight, NULL, 0.0, 0, 0, 0, ""); uiDefButF(block, NUM, B_DIFF, "", 150, yline,150,objHeight, &fss->particleInfSize, 0.0, 2.0, 10,0, "Amount of particle size scaling: 0=off (all same size), 1=full (range 0.2-2.0), >1=stronger."); @@ -2840,7 +2880,7 @@ static void object_panel_fluidsim(Object *ob) uiDefBut(block, TEX, B_FLUIDSIM_FORCEREDRAW,"", 20, yline, 280, objHeight, fss->surfdataPath, 0.0,79.0, 0, 0, "Enter fluid simulation bake directory/prefix to load particles from, same as for domain object."); uiBlockEndAlign(block); yline -= lineHeight; - + } // disabled for now... } else { diff --git a/source/blender/src/fluidsim.c b/source/blender/src/fluidsim.c index c555aea654e..2103d39a753 100644 --- a/source/blender/src/fluidsim.c +++ b/source/blender/src/fluidsim.c @@ -199,6 +199,7 @@ FluidsimSettings *fluidsimSettingsNew(struct Object *srcob) fss->bbSize[2] = 1.0; fluidsimGetAxisAlignedBB(srcob->data, srcob->obmat, fss->bbStart, fss->bbSize, &fss->meshBB); + // todo - reuse default init from elbeem! fss->typeFlags = 0; fss->domainNovecgen = 0; fss->volumeInitType = 1; // volume @@ -207,6 +208,7 @@ FluidsimSettings *fluidsimSettingsNew(struct Object *srcob) fss->generateTracers = 0; fss->generateParticles = 0.0; fss->surfaceSmoothing = 1.0; + fss->surfaceSubdivs = 1.0; fss->particleInfSize = 0.0; fss->particleInfAlpha = 0.0; @@ -915,6 +917,7 @@ void fluidsimBake(struct Object *ob) fsset.generateParticles = domainSettings->generateParticles; fsset.numTracerParticles = domainSettings->generateTracers; fsset.surfaceSmoothing = domainSettings->surfaceSmoothing; + fsset.surfaceSubdivs = domainSettings->surfaceSubdivs; fsset.farFieldSize = domainSettings->farFieldSize; strcpy( fsset.outputPath, targetFile); @@ -929,10 +932,10 @@ void fluidsimBake(struct Object *ob) fsset.runsimCallback = &runSimulationCallback; fsset.runsimUserData = &fsset; - if( (domainSettings->typeFlags&OB_FSBND_NOSLIP)) fsset.obstacleType = FLUIDSIM_OBSTACLE_NOSLIP; - 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; + if( (domainSettings->typeFlags&OB_FSBND_NOSLIP)) fsset.domainobsType = FLUIDSIM_OBSTACLE_NOSLIP; + else if((domainSettings->typeFlags&OB_FSBND_PARTSLIP)) fsset.domainobsType = FLUIDSIM_OBSTACLE_PARTSLIP; + else if((domainSettings->typeFlags&OB_FSBND_FREESLIP)) fsset.domainobsType = FLUIDSIM_OBSTACLE_FREESLIP; + fsset.domainobsPartslip = domainSettings->partSlipValue; fsset.generateVertexVectors = (domainSettings->domainNovecgen==0); // init blender trafo matrix @@ -962,7 +965,7 @@ void fluidsimBake(struct Object *ob) int *tris=NULL; int numVerts=0, numTris=0; int o = channelObjCount; - int deform = (obit->fluidsimSettings->domainNovecgen); + int deform = (obit->fluidsimSettings->domainNovecgen); // misused value elbeemMesh fsmesh; elbeemResetMesh( &fsmesh ); fsmesh.type = obit->fluidsimSettings->type;; @@ -996,6 +999,7 @@ void fluidsimBake(struct Object *ob) else if((obit->fluidsimSettings->typeFlags&OB_FSBND_FREESLIP)) fsmesh.obstacleType = FLUIDSIM_OBSTACLE_FREESLIP; fsmesh.obstaclePartslip = obit->fluidsimSettings->partSlipValue; fsmesh.volumeInitType = obit->fluidsimSettings->volumeInitType; + fsmesh.obstacleImpactFactor = obit->fluidsimSettings->surfaceSmoothing; // misused value // animated meshes if(deform) { @@ -1093,287 +1097,7 @@ void fluidsimBake(struct Object *ob) // -------------------------------------------------------------------------------------------- else { // write config file to be run with command line simulator - fileCfg = fopen(targetFile, "w"); - if(!fileCfg) { - snprintf(debugStrBuffer,256,"fluidsimBake::error - Unable to open file for writing '%s'\n", targetFile); - elbeemDebugOut(debugStrBuffer); - - pupmenu("Fluidsim Bake Error%t|Unable to output files... Aborted%x0"); - FS_FREE_CHANNELS; - return; - } - //ADD_CREATEDFILE(targetFile); - - fprintf(fileCfg, "# Blender ElBeem File , Source %s , Frame %d, to %s \n\n\n", G.sce, -1, targetFile ); - // file open -> valid settings -> store - strncpy(domainSettings->surfdataPath, newSurfdataPath, FILE_MAXDIR); - - /* output simulation settings */ - { - char *dtype[3] = { "no", "part", "free" }; - float pslip = domainSettings->partSlipValue; int bi=0; - char *simString = "\n" - "attribute \"simulation1\" { \n" - " solver = \"fsgr\"; \n" "\n" - " p_domainsize = " "%f" /* realsize */ "; \n" - " p_anistart = " "%f" /* aniStart*/ "; \n" - " p_normgstar = %f; \n" /* use gstar param? */ - " maxrefine = " "%d" /* maxRefine*/ "; \n" - " size = " "%d" /* gridSize*/ "; \n" - " surfacepreview = " "%d" /* previewSize*/ "; \n" - " dump_velocities = " "%d" /* vector dump */ "; \n" - " smoothsurface = %f; \n" /* smoothing */ - " domain_trafo = %f %f %f %f " /* remove blender object trafo */ - " %f %f %f %f " - " %f %f %f %f " - " %f %f %f %f ;\n" - " smoothnormals = %f; \n" - " geoinitid = 1; \n" "\n" - " isovalue = 0.4900; \n" - " isoweightmethod = 1; \n" "\n" ; - - fprintf(fileCfg, simString, - (double)domainSettings->realsize, (double)domainSettings->animStart, (double)domainSettings->gstar, - gridlevels, (int)domainSettings->resolutionxyz, (int)domainSettings->previewresxyz , - (int)(domainSettings->domainNovecgen==0), domainSettings->surfaceSmoothing, - invDomMat[0][0],invDomMat[1][0],invDomMat[2][0],invDomMat[3][0], - invDomMat[0][1],invDomMat[1][1],invDomMat[2][1],invDomMat[3][1], - invDomMat[0][2],invDomMat[1][2],invDomMat[2][2],invDomMat[3][2], - invDomMat[0][3],invDomMat[1][3],invDomMat[2][3],invDomMat[3][3] - ); - - if((domainSettings->typeFlags&OB_FSBND_NOSLIP)) bi=0; - else if((domainSettings->typeFlags&OB_FSBND_PARTSLIP)) bi=1; - else if((domainSettings->typeFlags&OB_FSBND_FREESLIP)) bi=2; - fprintf(fileCfg, " domainbound = %s; domainpartslip=%f; \n", dtype[bi], pslip); - - fprintf(fileCfg," # org aniframetime: %f \n", aniFrameTime); - fluidsimPrintChannel(fileCfg, channelDomainTime,allchannelSize,"p_aniframetime",CHANNEL_FLOAT); - 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, " particles = %d; \n", domainSettings->generateTracers); // debug test - fprintf(fileCfg, "\n} \n" ); - } - - - fprintf(fileCfg, "raytracing {\n"); - - /* output picture settings for preview renders */ - { - char *rayString = "\n" - " anistart= 0; \n" - " aniframes= " "%d" /*1 frameEnd-frameStart+0*/ "; #cfgset \n" - " frameSkip= false; \n" - " filename= \"" "%s" /* rayPicFilename*/ "\"; #cfgset \n" - " aspect 1.0; \n" - " resolution " "%d %d" /*2,3 blendResx,blendResy*/ "; #cfgset \n" - " antialias 1; \n" - " ambientlight (1, 1, 1); \n" - " maxRayDepth 6; \n" - " treeMaxDepth 25; \n" - " treeMaxTriangles 8; \n" - " background (0.08, 0.08, 0.20); \n" - " eyepoint= (" "%f %f %f"/*4,5,6 eyep*/ "); #cfgset \n" - " lookat= (" "%f %f %f"/*7,8,9 lookatp*/ "); #cfgset \n" - " upvec= (0 0 1); \n" - " fovy= " "%f" /*blendFov*/ "; #cfgset \n" - //" blenderattr= \"btrafoattr\"; \n" - "\n\n"; - - char *lightString = "\n" - " light { \n" - " type= omni; \n" - " active= 1; \n" - " color= (1.0, 1.0, 1.0); \n" - " position= (" "%f %f %f"/*1,2,3 eyep*/ "); #cfgset \n" - " castShadows= 1; \n" - " } \n\n" ; - - struct Object *cam = G.scene->camera; - float eyex=2.0, eyey=2.0, eyez=2.0; - int resx = 200, resy=200; - float lookatx=0.0, lookaty=0.0, lookatz=0.0; - float fov = 45.0; - - strcpy(targetFile, targetDir); - strcat(targetFile, suffixSurface); - resx = G.scene->r.xsch; - resy = G.scene->r.ysch; - if((cam) && (cam->type == OB_CAMERA)) { - Camera *camdata= G.scene->camera->data; - double lens = camdata->lens; - double imgRatio = (double)resx/(double)resy; - fov = 360.0 * atan(16.0*imgRatio/lens) / M_PI; - //R.near= camdata->clipsta; R.far= camdata->clipend; - - eyex = cam->loc[0]; - eyey = cam->loc[1]; - eyez = cam->loc[2]; - // TODO - place lookat in middle of domain? - } - - fprintf(fileCfg, rayString, - (noFrames+0), targetFile, resx,resy, - eyex, eyey, eyez , - lookatx, lookaty, lookatz, - fov - ); - fprintf(fileCfg, lightString, - eyex, eyey, eyez ); - } - - - /* output fluid domain */ - { - char * domainString = "\n" - " geometry { \n" - " type= fluidlbm; \n" - " name = \"" "%s" /*name*/ "\"; #cfgset \n" - " visible= 1; \n" - " attributes= \"simulation1\"; \n" - //" define { material_surf = \"fluidblue\"; } \n" - " start= " "%f %f %f" /*bbstart*/ "; #cfgset \n" - " end = " "%f %f %f" /*bbend */ "; #cfgset \n" - " } \n" - "\n"; - fprintf(fileCfg, domainString, - fsDomain->id.name, - bbStart[0], bbStart[1], bbStart[2], - bbStart[0]+bbSize[0], bbStart[1]+bbSize[1], bbStart[2]+bbSize[2] - ); - } - - - /* setup geometry */ - { - char *objectStringStart = - " geometry { \n" - " type= objmodel; \n" - " name = \"" "%s" /* name */ "\"; #cfgset \n" - " visible= 1; \n" // DEBUG , also obs invisible? - " define { \n" ; - char *outflowString = - " geoinittype= \"" "%s" /* type */ "\"; #cfgset \n" - " filename= \"" "%s" /* data filename */ "\"; #cfgset \n" ; - char *obstacleString = - " geoinittype= \"" "%s" /* type */ "\"; #cfgset \n" - " geoinit_partslip = \"" "%f" /* partslip */ "\"; #cfgset \n" - " geoinit_volumeinit = \"" "%d" /* volumeinit */ "\"; #cfgset \n" - " filename= \"" "%s" /* data filename */ "\"; #cfgset \n" ; - char *fluidString = - " geoinittype= \"" "%s" /* type */ "\"; \n" - " geoinit_volumeinit = \"" "%d" /* volumeinit */ "\"; #cfgset \n" - " filename= \"" "%s" /* data filename */ "\"; #cfgset \n" ; - char *inflowString = - " geoinittype= \"" "%s" /* type */ "\"; \n" - " geoinit_localinivel = " "%d" /* local coords */ "; #cfgset \n" - " filename= \"" "%s" /* data filename */ "\"; #cfgset \n" ; - char *objectStringEnd = - " geoinit_intersect = 1; \n" /* always use accurate init here */ - " geoinitid= 1; \n" - " } \n" - " } \n" - "\n" ; - char fnameObjdat[FILE_MAXFILE]; - - channelObjCount = 0; - for(obit= G.main->object.first; obit; obit= obit->id.next) { - if( (obit->fluidsimFlag & OB_FLUIDSIM_ENABLE) && - (obit->type==OB_MESH) && // if has to match 3 places! // CHECKMATCH - (obit->fluidsimSettings->type != OB_FLUIDSIM_DOMAIN) && - (obit->fluidsimSettings->type != OB_FLUIDSIM_PARTICLE) - ) { - int deform = (obit->fluidsimSettings->domainNovecgen); - - fluidsimGetGeometryObjFilename(obit, fnameObjdat); - strcpy(targetFile, targetDir); - strcat(targetFile, fnameObjdat); - fprintf(fileCfg, objectStringStart, obit->id.name ); // abs path - // object type params - if(obit->fluidsimSettings->type == OB_FLUIDSIM_FLUID) { - fprintf(fileCfg, fluidString, "fluid", - (int)obit->fluidsimSettings->volumeInitType, targetFile ); - } - if(obit->fluidsimSettings->type == OB_FLUIDSIM_INFLOW) { - int locc = ((obit->fluidsimSettings->typeFlags&OB_FSINFLOW_LOCALCOORD)?1:0); - fprintf(fileCfg, inflowString, "inflow" ,locc - , targetFile ); - } - if(obit->fluidsimSettings->type == OB_FLUIDSIM_OBSTACLE) { - char *btype[3] = { "bnd_no", "bnd_part", "bnd_free" }; - float pslip = obit->fluidsimSettings->partSlipValue; int bi=0; - if((obit->fluidsimSettings->typeFlags&OB_FSBND_NOSLIP)) bi=0; - else if((obit->fluidsimSettings->typeFlags&OB_FSBND_PARTSLIP)) bi=1; - else if((obit->fluidsimSettings->typeFlags&OB_FSBND_FREESLIP)) bi=2; - fprintf(fileCfg, obstacleString, btype[bi], pslip, - (int)obit->fluidsimSettings->volumeInitType, targetFile); // abs path - } - if(obit->fluidsimSettings->type == OB_FLUIDSIM_OUTFLOW) { - fprintf(fileCfg, outflowString, "outflow" , targetFile); // abs path - } - - if(!deform) { - fluidsimPrintChannel(fileCfg, channelObjMove[channelObjCount][0],allchannelSize, "translation", CHANNEL_VEC); - fluidsimPrintChannel(fileCfg, channelObjMove[channelObjCount][1],allchannelSize, "rotation" , CHANNEL_VEC); - fluidsimPrintChannel(fileCfg, channelObjMove[channelObjCount][2],allchannelSize, "scale" , CHANNEL_VEC); - } - fluidsimPrintChannel(fileCfg, channelObjActive[channelObjCount] ,allchannelSize, "geoactive" , CHANNEL_FLOAT); - if( (obit->fluidsimSettings->type == OB_FLUIDSIM_FLUID) || - (obit->fluidsimSettings->type == OB_FLUIDSIM_INFLOW) ) { - fluidsimPrintChannel(fileCfg, channelObjInivel[channelObjCount],allchannelSize,"initial_velocity" ,CHANNEL_VEC); - } - channelObjCount++; - - fprintf(fileCfg, objectStringEnd ); // abs path - - // check shape key animation - //fprintf(stderr,"\n%d %d\n\n",(int)obit->parent,obit->partype); // DEBUG - if(deform) { - int frame; - // use global coordinates for deforming/parented objects - writeBobjgz(targetFile, obit, 1,0,0.); - //for(int frame=0; frame<=G.scene->r.efra; frame++) { - for(frame=0; frame<=allchannelSize; frame++) { - G.scene->r.cfra = frame; - scene_update_for_newframe(G.scene, G.scene->lay); - writeBobjgz(targetFile, obit, 1,1, timeAtFrame[frame] ); // only append! - - //if(shapekey) snprintf(debugStrBuffer,256,"Shape frames: %d/%d, shapeKeys:%d",frame,allchannelSize,BLI_countlist(&shapekey->block)); - //else snprintf(debugStrBuffer,256,"Deform frames: %d/%d",frame,allchannelSize); - //elbeemDebugOut(debugStrBuffer); - } - G.scene->r.cfra = startFrame; - scene_update_for_newframe(G.scene, G.scene->lay); - } else { - // use normal trafos & non animated mesh - writeBobjgz(targetFile, obit, 0,0,0.); - } - - } - } - } - - /* fluid material */ - fprintf(fileCfg, - " material { \n" - " type= phong; \n" - " name= \"fluidblue\"; \n" - " diffuse= 0.3 0.5 0.9; \n" - " ambient= 0.1 0.1 0.1; \n" - " specular= 0.2 10.0; \n" - " } \n" ); - - fprintf(fileCfg, "} // end raytracing\n"); - fclose(fileCfg); - - strcpy(targetFile, targetDir); - strcat(targetFile, suffixConfig); - snprintf(debugStrBuffer,256,"fluidsimBake::msg: Wrote %s\n", targetFile); - elbeemDebugOut(debugStrBuffer); - - pupmenu("Fluidsim Bake Message%t|Config files exported successfully!%x0"); + pupmenu("Fluidsim Bake Message%t|Config file export not supported.%x0"); } // config file export done! // --------------------------------------------------------------------------------------------