archive/blender-archive
Archived
1
1
Fork 0
Code Pull Requests Packages Activity

- bugfixes

Browse Source 
#4742 exported normals are now correct
  #4821 & 4956 for complex movements in/outflows can now also
  use the animated mesh option
- new features
  * isosurface subdivision: directly
    creates a finer surface mesh from the simulation data.
    this increases simulation time and harddisk usage, though, so
    be careful - usually values of 2-4 should be enough.
  * fluidsim particles: extended model for particle
    simulation and generation. When isosurface subdivision is enabled,
    the particles are now included in the surface generation,
    giving a better impression of a single connected surface.
    Note - the particles are only included in the final surface
    mesh, so the preview surface shows none of the particle
    effects.
  * particle loading: different types of particles can now be selected for
    display: drops, floats and tracers. This is a bit obsolete
    due to the extensions mentioned above, but might still be useful.
    Floats are just particles floating on the fluid surface, could
    be used for e.g. foam.
  * moving objects impact factor: this is another tweaking option,
    as the handling of moving objects is still not conserving
    mass. setting this to zero simply deletes the fluid, 1 is
    the default, while larger values cause a stronger
    impact. For tweaking the simulation: if fluid disappears, try
    increasing this value, and if too much is appearing reduce it.
    You can even use negative values for some strange results :)
- more code cleanup, e.g. removed config file writing in fluidsim.c,
  added additional safety checks for particles & fluidsim domains (these
  currently dont work together). I also removed the "build particles"
  debug message in effects.c (seemed to be unnecessary?).

Some more info on the new features:
Here are two test animations showing the difference between
using the particle generation with isosurface subdivision.
This is how it would look with the old solver version:
http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid6_fl6manc4_1noparts.mpg
and this with the new one:
http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid6_fl6manc4_2wparts.mpg
Both simulations use a resolution of 64, however, the version with particles
takes significantly longer (almost twice as long).
The .blend file for a similar setup can be found here:
http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid6_testmanc4.blend
(Minor Tips for this file: dont enable subdivions of characters until rendering,
thus leave off for simulation, as it uses the rendering settings! For making
nice pictures switch on subdivion, and OSA.)

And here's a picture of old vs. new (for webpage or so):
http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid6_manc4compare.png
This commit is contained in:
Nils Thuerey
2006-11-05 16:30:29 +00:00
parent 64b9cda68e
commit 3bea663ffa
48 changed files with 2601 additions and 2345 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
intern/elbeem/COPYING
View File

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

14
intern/elbeem/extern/elbeem.h vendored
View File

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

556
intern/elbeem/intern/attributes.cpp
View File

@@ -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 <sstream>
//! 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!? ("<<mChannelInited<<" vs "<<elemSize<<")...");
return false;
}
}
if((mValue.size() % (elemSize+1)) != 0) {
errMsg("Attribute::initChannel","Invalid # elements in Attribute...");
return false;
}
int numElems = mValue.size()/(elemSize+1);
vector<string> newvalue;
for(int i=0; i<numElems; i++) {
//errMsg("ATTR"," i"<<i<<" "<<mName); // debug
vector<string> elem(elemSize);
for(int j=0; j<elemSize; j++) {
//errMsg("ATTR"," j"<<j<<" "<<mValue[i*(elemSize+1)+j] ); // debug
elem[j] = mValue[i*(elemSize+1)+j];
}
mChannel.push_back(elem);
// use first value as default
if(i==0) newvalue = elem;
double t = 0.0; // similar to getAsFloat
const char *str = mValue[i*(elemSize+1)+elemSize].c_str();
char *endptr;
t = strtod(str, &endptr);
if((str!=endptr) && (*endptr != '\0')) return false;
mTimes.push_back(t);
//errMsg("ATTR"," t"<<t<<" "); // debug
}
for(int i=0; i<numElems-1; i++) {
if(mTimes[i]>mTimes[i+1]) {
errMsg("Attribute::initChannel","Invalid time at entry "<<i<<" setting to "<<mTimes[i]);
mTimes[i+1] = mTimes[i];
}
}
// dont change until done with parsing, and everythings ok
mValue = newvalue;
mChannelInited = elemSize;
if(DEBUG_CHANNELS) print();
return true;
return false;
}
// get value as string
string Attribute::getAsString(bool debug)
{
if(mIsChannel && (!debug)) {
errMsg("Attribute::getAsString", "Attribute \"" << mName << "\" used as string is a channel! Not allowed...");
print();
return string("");
}
if(mValue.size()!=1) {
// for directories etc. , this might be valid! cutoff "..." first
string comp = getCompleteString();
if(comp.size()<2) return string("");
return comp.substr(1, comp.size()-2);
}
return mValue[0];
string Attribute::getAsString(bool debug) {
return string("");
}
// get value as integer value
int Attribute::getAsInt()
{
bool success = true;
int ret = 0;
if(!initChannel(1)) success=false;
if(success) {
if(mValue.size()!=1) success = false;
else {
const char *str = mValue[0].c_str();
char *endptr;
ret = strtol(str, &endptr, 10);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
}
}
if(!success) {
errMsg("Attribute::getAsInt", "Attribute \"" << mName << "\" used as int has invalid value '"<< getCompleteString() <<"' ");
errMsg("Attribute::getAsInt", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
errMsg("Attribute::getAsInt", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
#if ELBEEM_PLUGIN!=1
setElbeemState( -4 ); // parse error
#endif
return 0;
}
return ret;
int Attribute::getAsInt() {
return 0;
}
// get value as integer value
bool Attribute::getAsBool()
{
int val = getAsInt();
if(val==0) return false;
else return true;
bool Attribute::getAsBool() {
return false;
}
// get value as double value
double Attribute::getAsFloat()
{
bool success = true;
double ret = 0.0;
if(!initChannel(1)) success=false;
if(success) {
if(mValue.size()!=1) success = false;
else {
const char *str = mValue[0].c_str();
char *endptr;
ret = strtod(str, &endptr);
if((str!=endptr) && (*endptr != '\0')) success = false;
}
}
if(!success) {
print();
errMsg("Attribute::getAsFloat", "Attribute \"" << mName << "\" used as double has invalid value '"<< getCompleteString() <<"' ");
errMsg("Attribute::getAsFloat", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
errMsg("Attribute::getAsFloat", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
#if ELBEEM_PLUGIN!=1
setElbeemState( -4 ); // parse error
#endif
return 0.0;
}
return ret;
double Attribute::getAsFloat() {
return 0.;
}
// get value as 3d vector
ntlVec3d Attribute::getAsVec3d()
{
bool success = true;
ntlVec3d ret(0.0);
if(!initChannel(3)) success=false;
if(success) {
if(mValue.size()==1) {
const char *str = mValue[0].c_str();
char *endptr;
double rval = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
if(success) ret = ntlVec3d( rval );
} else if(mValue.size()==3) {
char *endptr;
const char *str = NULL;
str = mValue[0].c_str();
double rval1 = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
str = mValue[1].c_str();
double rval2 = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
str = mValue[2].c_str();
double rval3 = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
if(success) ret = ntlVec3d( rval1, rval2, rval3 );
} else {
success = false;
}
}
if(!success) {
errMsg("Attribute::getAsVec3d", "Attribute \"" << mName << "\" used as Vec3d has invalid value '"<< getCompleteString() <<"' ");
errMsg("Attribute::getAsVec3d", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
errMsg("Attribute::getAsVec3d", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
#if ELBEEM_PLUGIN!=1
setElbeemState( -4 ); // parse error
#endif
return ntlVec3d(0.0);
}
return ret;
ntlVec3d Attribute::getAsVec3d() {
return ntlVec3d(0.);
}
// get value as 4x4 matrix
void Attribute::getAsMat4Gfx(ntlMat4Gfx *mat)
{
bool success = true;
ntlMat4Gfx ret(0.0);
char *endptr;
if(mValue.size()==1) {
const char *str = mValue[0].c_str();
double rval = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) )success = false;
if(success) {
ret = ntlMat4Gfx( 0.0 );
ret.value[0][0] = rval;
ret.value[1][1] = rval;
ret.value[2][2] = rval;
ret.value[3][3] = 1.0;
}
} else if(mValue.size()==9) {
// 3x3
for(int i=0; i<3;i++) {
for(int j=0; j<3;j++) {
const char *str = mValue[i*3+j].c_str();
ret.value[i][j] = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) ) success = false;
}
}
} else if(mValue.size()==16) {
// 4x4
for(int i=0; i<4;i++) {
for(int j=0; j<4;j++) {
const char *str = mValue[i*4+j].c_str();
ret.value[i][j] = strtod(str, &endptr);
if( (str==endptr) ||
((str!=endptr) && (*endptr != '\0')) ) success = false;
}
}
} else {
success = false;
}
if(!success) {
errMsg("Attribute::getAsMat4Gfx", "Attribute \"" << mName << "\" used as Mat4x4 has invalid value '"<< getCompleteString() <<"' ");
errMsg("Attribute::getAsMat4Gfx", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
errMsg("Attribute::getAsMat4Gfx", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
#if ELBEEM_PLUGIN!=1
setElbeemState( -4 ); // parse error
#endif
*mat = ntlMat4Gfx(0.0);
return;
}
*mat = ret;
void Attribute::getAsMat4Gfx(ntlMat4Gfx *mat) {
}
// get the concatenated string of all value string
string Attribute::getCompleteString()
{
string ret;
for(size_t i=0;i<mValue.size();i++) {
ret += mValue[i];
if(i<mValue.size()-1) ret += " ";
}
return ret;
string Attribute::getCompleteString() {
return string("");
}
@@ -292,266 +46,80 @@ string Attribute::getCompleteString()
* channel returns
*****************************************************************************/
//! get channel as double value
AnimChannel<double> Attribute::getChannelFloat() {
vector<double> timeVec;
vector<double> valVec;
if((!initChannel(1)) || (!mIsChannel)) {
timeVec.push_back( 0.0 );
valVec.push_back( getAsFloat() );
} else {
for(size_t i=0; i<mChannel.size(); i++) {
mValue = mChannel[i];
double val = getAsFloat();
timeVec.push_back( mTimes[i] );
valVec.push_back( val );
}}
return AnimChannel<double>(valVec,timeVec);
return AnimChannel<double>();
}
//! get channel as integer value
AnimChannel<int> Attribute::getChannelInt() {
vector<double> timeVec;
vector<int> valVec;
if((!initChannel(1)) || (!mIsChannel)) {
timeVec.push_back( 0.0 );
valVec.push_back( getAsInt() );
} else {
for(size_t i=0; i<mChannel.size(); i++) {
mValue = mChannel[i];
int val = getAsInt();
timeVec.push_back( mTimes[i] );
valVec.push_back( val );
}}
return AnimChannel<int>(valVec,timeVec);
return AnimChannel<int>();
}
//! get channel as integer value
AnimChannel<ntlVec3d> Attribute::getChannelVec3d() {
vector<double> timeVec;
vector<ntlVec3d> valVec;
if((!initChannel(3)) || (!mIsChannel)) {
timeVec.push_back( 0.0 );
valVec.push_back( getAsVec3d() );
} else {
for(size_t i=0; i<mChannel.size(); i++) {
mValue = mChannel[i];
ntlVec3d val = getAsVec3d();
timeVec.push_back( mTimes[i] );
valVec.push_back( val );
}}
return AnimChannel<ntlVec3d>(valVec,timeVec);
return AnimChannel<ntlVec3d>();
}
//! get channel as float vector set
AnimChannel<ntlSetVec3f>
Attribute::getChannelSetVec3f() {
vector<double> 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<mChannel.size(); i++) {
mValue = mChannel[i];
ntlVec3f val = vec2F( getAsVec3d() );
timeVec.push_back( mTimes[i] );
setv.mVerts.push_back( val );
}
valVec.push_back( setv );
valVec.push_back( setv );
}
return AnimChannel<ntlSetVec3f>(valVec,timeVec);
return AnimChannel<ntlSetVec3f>();
}
/******************************************************************************
* check if there were unknown params
*****************************************************************************/
bool AttributeList::checkUnusedParams()
{
bool found = false;
for(map<string, Attribute*>::iterator i=mAttrs.begin();
i != mAttrs.end(); i++) {
if((*i).second) {
if(!(*i).second->getUsed()) {
errMsg("AttributeList::checkUnusedParams", "List "<<mName<<" has unknown parameter '"<<(*i).first<<"' = '"<< mAttrs[(*i).first]->getAsString(true) <<"' ");
found = true;
}
}
}
return found;
bool AttributeList::checkUnusedParams() {
return false;
}
//! set all params to used, for invisible objects
void AttributeList::setAllUsed() {
for(map<string, Attribute*>::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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
return defaultValue;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
return defaultValue;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
return defaultValue;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
return defaultValue;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
return defaultValue;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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 '"<<name<<"' for "<< source <<" not set! ", SIMWORLD_INITERROR); }
*mat = defaultValue;
return;
}
if(DEBUG_ATTRIBUTES==1) { debugOut( source << " Var '"<< target <<"' set to '"<< find(name)->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<int> AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<int>(defaultValue); }
AnimChannel<int> ret = find(name)->getChannelInt();
find(name)->setUsed(true);
channelSimplifyi(ret);
return ret;
return AnimChannel<int>(defaultValue);
}
AnimChannel<double> AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<double>(defaultValue); }
AnimChannel<double> ret = find(name)->getChannelFloat();
find(name)->setUsed(true);
channelSimplifyd(ret);
return ret;
return AnimChannel<double>(defaultValue);
}
AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<ntlVec3d>(defaultValue); }
AnimChannel<ntlVec3d> ret = find(name)->getChannelVec3d();
find(name)->setUsed(true);
channelSimplifyVd(ret);
return ret;
return AnimChannel<ntlVec3d>(defaultValue);
}
AnimChannel<ntlSetVec3f> AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<ntlSetVec3f>(defaultValue); }
AnimChannel<ntlSetVec3f> ret = find(name)->getChannelSetVec3f();
find(name)->setUsed(true);
//channelSimplifyVf(ret);
return ret;
return AnimChannel<ntlSetVec3f>(defaultValue);
}
AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<float>(defaultValue); }
AnimChannel<double> convert = find(name)->getChannelFloat();
find(name)->setUsed(true);
channelSimplifyd(convert);
// convert to float vals
vector<float> vals;
for(size_t i=0; i<convert.accessValues().size(); i++) {
vals.push_back( (float)(convert.accessValues()[i]) );
}
vector<double> times = convert.accessTimes();
AnimChannel<float> ret(vals, times);
return ret;
return AnimChannel<float>(defaultValue);
}
AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) {
if(!exists(name)) {
if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<" not set! ", SIMWORLD_INITERROR); }
return AnimChannel<ntlVec3f>(defaultValue); }
AnimChannel<ntlVec3d> convert = find(name)->getChannelVec3d();
// convert to float
vector<ntlVec3f> vals;
for(size_t i=0; i<convert.accessValues().size(); i++) {
vals.push_back( vec2F(convert.accessValues()[i]) );
}
vector<double> times = convert.accessTimes();
AnimChannel<ntlVec3f> ret(vals, times);
find(name)->setUsed(true);
channelSimplifyVf(ret);
return ret;
return AnimChannel<ntlVec3f>(defaultValue);
}
/******************************************************************************
* destructor
*****************************************************************************/
AttributeList::~AttributeList() {
for(map<string, Attribute*>::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;i<mValue.size();i++) {
ostr <<"'"<< mValue[i]<<"' ";
}
if(mIsChannel) {
ostr << " CHANNEL: ";
if(mChannelInited>0) {
for(size_t i=0;i<mChannel.size();i++) {
for(size_t j=0;j<mChannel[i].size();j++) {
ostr <<"'"<< mChannel[i][j]<<"' ";
}
ostr << "@"<<mTimes[i]<<"; ";
}
} else {
ostr <<" -nyi- ";
}
}
ostr <<" (at line "<<mLine<<") "; //<< std::endl;
debugOut( ostr.str(), 10);
void Attribute::print() {
}
//! debug function, prints all attribs
void AttributeList::print()
{
debugOut("Attribute "<<mName<<" values:", 10);
for(map<string, Attribute*>::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<string, Attribute*>::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<SCALAR> &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.accessValues().size(); i++) {
nchannel[i*2 + 0] = (float)channel.accessValues()[i];
@@ -687,7 +216,6 @@ static bool channelSimplifyScalarT(AnimChannel<SCALAR> &channel) {
times.push_back( (double)(nchannel[i*2 + 1]) );
}
channel = AnimChannel<SCALAR>(vals, times);
if(DEBUG_CHANNELS) errMsg("channelSimplifyf","C" << channel.printChannel() );
}
delete [] nchannel;
return ret;
@@ -700,7 +228,6 @@ static bool channelSimplifyVecT(AnimChannel<VEC> &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.accessValues().size(); i++) {
nchannel[i*4 + 0] = (float)channel.accessValues()[i][0];
@@ -717,7 +244,6 @@ static bool channelSimplifyVecT(AnimChannel<VEC> &channel) {
times.push_back( (double)(nchannel[i*4 + 3]) );
}
channel = AnimChannel<VEC>(vals, times);
if(DEBUG_CHANNELS) errMsg("channelSimplifyf","C" << channel.printChannel() );
}
delete [] nchannel;
return ret;
@@ -742,13 +268,7 @@ string AnimChannel<Scalar>::printChannel() {
return ostr.str();
} // */
//! debug function, prints to stdout if DEBUG_CHANNELS flag is enabled, used in constructors
template<class Scalar>
void AnimChannel<Scalar>::debugPrintChannel() {
if(DEBUG_CHANNELS) { errMsg("channelCons"," : " << this->printChannel() ); }
}
// is now in header file: debugPrintChannel()
// hack to force instantiation
void __forceAnimChannelInstantiation() {
AnimChannel< float > tmp1;

192
intern/elbeem/intern/attributes.h
View File

@@ -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"<<t<<" ");
if(!mInited) { Scalar null; null=(Scalar)0.0; return null; }
if(t<=mTimes[0]) { return mValue[0]; }
@@ -90,10 +90,10 @@ class AnimChannel
//! debug function, prints to stdout if DEBUG_CHANNELS flag is enabled, used in constructors
void debugPrintChannel();
//! valid init?
bool isInited() { return mInited; }
bool isInited() const { return mInited; }
//! get number of entries (value and time sizes have to be equal)
int getSize() { return mValue.size(); };
int getSize() const { return mValue.size(); };
//! raw access of value vector
vector<Scalar> &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<string> &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<gfxReal> *mat);
//! get channel as integer value
AnimChannel<int> getChannelInt();
//! get channel as double value
AnimChannel<double> getChannelFloat();
//! get channel as double vector
AnimChannel<ntlVec3d> getChannelVec3d();
//! get channel as float vector set
AnimChannel<ntlSetVec3f> 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<string> 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<ntlVec3f> mVerts;
};
// warning: DEPRECATED - replaced by elbeem API
class Attribute
{
public:
Attribute(string mn, vector<string> &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<gfxReal> *mat);
AnimChannel<int> getChannelInt();
AnimChannel<double> getChannelFloat();
AnimChannel<ntlVec3d> getChannelVec3d();
AnimChannel<ntlSetVec3f> 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<string> &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 '"<<name<<"' , not found...",SIMWORLD_INITERROR );
// just create a new empty one (warning: small memory leak!), and exit as soon as possible
vector<string> empty;
return new Attribute(name,empty, -1, 0);
}
return mAttrs[name];
}
//! set all params to used, for invisible objects
void addAttr(string name, vector<string> &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<gfxReal> defaultValue, string source,string target, bool needed, ntlMatrix4x4<gfxReal> *mat);
//! read attributes channels (attribute should be inited before)
AnimChannel<int> readChannelInt( string name, int defaultValue=0, string source=string("src"), string target=string("dst"), bool needed=false );
AnimChannel<double> readChannelFloat( string name, double defaultValue=0, string source=string("src"), string target=string("dst"), bool needed=false );
AnimChannel<ntlVec3d> readChannelVec3d( string name, ntlVec3d defaultValue=ntlVec3d(0.), string source=string("src"), string target=string("dst"), bool needed=false );
AnimChannel<ntlSetVec3f> readChannelSetVec3f(string name, ntlSetVec3f defaultValue=ntlSetVec3f(0.), string source=string("src"), string target=string("dst"), bool needed=false );
// channels with conversion
AnimChannel<ntlVec3f> readChannelVec3f( string name, ntlVec3f defaultValue=ntlVec3f(0.), string source=string("src"), string target=string("dst"), bool needed=false );
AnimChannel<float> 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<string, Attribute*> mAttrs;
};
ntlVec3f channelFindMaxVf (AnimChannel<ntlVec3f> channel);
@@ -297,6 +203,14 @@ bool channelSimplifyi (AnimChannel<int > &channel);
bool channelSimplifyf (AnimChannel<float > &channel);
bool channelSimplifyd (AnimChannel<double > &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<class Scalar>
void AnimChannel<Scalar>::debugPrintChannel() { }
#define NTL_ATTRIBUTES_H
#endif

10
intern/elbeem/intern/elbeem.cpp
View File

@@ -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:"<<gDebugLevel<<" ...\n", 2);
@@ -150,6 +152,7 @@ void elbeemResetMesh(elbeemMesh *mesh) {
mesh->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->volumeInitType<VOLUMEINIT_VOLUME)||(mesh->volumeInitType>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]) );

14
intern/elbeem/intern/elbeem.h
View File

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

730
intern/elbeem/intern/isosurface.cpp
View File

@@ -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 <algorithm>
#include <stdio.h>
// sirdude fix for solaris
#if !defined(linux) && (defined (__sparc) || defined (__sparc__))
#include <ieeefp.h>
@@ -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;i<nodes;i++) { mpData[i] = 0.0; }
// allocate edge arrays (last slices are never used...)
mpEdgeVerticesX = new int[nodes];
mpEdgeVerticesY = new int[nodes];
mpEdgeVerticesZ = new int[nodes];
for(int i=0;i<nodes;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
int initsize = -1;
if(mUseFullEdgeArrays) {
mEdgeArSize = nodes;
mpEdgeVerticesX = new int[nodes];
mpEdgeVerticesY = new int[nodes];
mpEdgeVerticesZ = new int[nodes];
initsize = 3*nodes;
} else {
int sliceNodes = 2*mSizex*mSizey*mSubdivs*mSubdivs;
mEdgeArSize = sliceNodes;
mpEdgeVerticesX = new int[sliceNodes];
mpEdgeVerticesY = new int[sliceNodes];
mpEdgeVerticesZ = new int[sliceNodes];
initsize = 3*sliceNodes;
}
for(int i=0;i<mEdgeArSize;i++) { mpEdgeVerticesX[i] = mpEdgeVerticesY[i] = mpEdgeVerticesZ[i] = -1; }
// WARNING - make sure this is consistent with calculateMemreqEstimate
// marching cubes are ready
mInitDone = true;
debMsgStd("IsoSurface::initializeIsosurface",DM_MSG,"Inited, edgenodes:"<<initsize<<" subdivs:"<<mSubdivs<<" fulledg:"<<mUseFullEdgeArrays , 10);
}
@@ -106,6 +124,10 @@ IsoSurface::~IsoSurface( void )
//static inline getSubdivData(IsoSurface* iso,int ai,aj,ak, si,sj) {
//int srci =
//}
/******************************************************************************
* triangulate the scalar field given by pointer
@@ -132,7 +154,7 @@ void IsoSurface::triangulate( void )
mPoints.clear();
// reset edge vertices
for(int i=0;i<(mSizex*mSizey*mSizez);i++) {
for(int i=0;i<mEdgeArSize;i++) {
mpEdgeVerticesX[i] = -1;
mpEdgeVerticesY[i] = -1;
mpEdgeVerticesZ[i] = -1;
@@ -160,150 +182,473 @@ void IsoSurface::triangulate( void )
const int coAdd=2;
// let the cubes march
pz = mStart[2]-gsz*0.5;
for(int k=1;k<(mSizez-2);k++) {
pz += gsz;
py = mStart[1]-gsy*0.5;
for(int j=1;j<(mSizey-2);j++) {
py += gsy;
px = mStart[0]-gsx*0.5;
for(int i=1;i<(mSizex-2);i++) {
px += gsx;
int baseIn = ISOLEVEL_INDEX( i+0, j+0, k+0);
if(mSubdivs<=1) {
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);
pz = mStart[2]-gsz*0.5;
for(int k=1;k<(mSizez-2);k++) {
pz += gsz;
py = mStart[1]-gsy*0.5;
for(int j=1;j<(mSizey-2);j++) {
py += gsy;
px = mStart[0]-gsx*0.5;
for(int i=1;i<(mSizex-2);i++) {
px += gsx;
/*int bndskip = 0; // BNDOFFT
for(int s=0; s<8; s++) if(value[s]==-76.) bndskip++;
if(bndskip>0) 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<<e)) {
// is the vertex already calculated?
if(*eVert[ e ] < 0) {
// interpolate edge
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1
const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2
const float valp1 = value[ e1 ]; // scalar field val 1
const float valp2 = value[ e2 ]; // scalar field val 2
float mu;
if(valp1 < valp2) {
mu = 1.0 - 1.0*(valp1 + valp2 - mIsoValue);
} else {
mu = 0.0 + 1.0*(valp1 + valp2 - mIsoValue);
}
//float isov2 = mIsoValue;
//isov2 = (valp1+valp2)*0.5;
//mu = (isov2 - valp1) / (valp2 - valp1);
//mu = (isov2) / (valp2 - valp1);
mu = (mIsoValue - valp1) / (valp2 - valp1);
// init isolevel vertex
ilv.v = p1 + (p2-p1)*mu;
ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * ( mu) ;
mPoints.push_back( ilv );
triIndices[e] = (mPoints.size()-1);
// store vertex
*eVert[ e ] = triIndices[e];
} else {
// retrieve from vert array
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
if( (i<coAdd+mCutoff) ||
(j<coAdd+mCutoff) ||
((mCutoff>0) && (k<coAdd)) ||// bottom layer
(i>mSizex-2-coAdd-mCutoff) ||
(j>mSizey-2-coAdd-mCutoff) ) {
if(mCutArray) {
if(k < mCutArray[j*this->mSizex+i]) continue;
} else {
// 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 "<<mIndices.size() <<" "<< triIndices[ mcTriTable[cubeIndex][e+0] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+1] ]<<" "<< triIndices[ mcTriTable[cubeIndex][e+2] ] );
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
} // k
// where to look up if this point already exists
int edgek = 0;
if(mUseFullEdgeArrays) edgek=k;
const int baseIn = ISOLEVEL_INDEX( i+0, j+0, edgek+0);
eVert[ 0] = &mpEdgeVerticesX[ baseIn ];
eVert[ 1] = &mpEdgeVerticesY[ baseIn + 1 ];
eVert[ 2] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+0) ];
eVert[ 3] = &mpEdgeVerticesY[ baseIn ];
// precalculate normals using an approximation of the scalar field gradient
for(int ni=0;ni<(int)mPoints.size();ni++) {
// use triangle normals?, this seems better for File-IsoSurf
normalize( mPoints[ni].n );
eVert[ 4] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 5] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+1, j+0, edgek+1) ];
eVert[ 6] = &mpEdgeVerticesX[ ISOLEVEL_INDEX( i+0, j+1, edgek+1) ];
eVert[ 7] = &mpEdgeVerticesY[ ISOLEVEL_INDEX( i+0, j+0, edgek+1) ];
eVert[ 8] = &mpEdgeVerticesZ[ baseIn ];
eVert[ 9] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+0, edgek+0) ];
eVert[10] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+1, j+1, edgek+0) ];
eVert[11] = &mpEdgeVerticesZ[ ISOLEVEL_INDEX( i+0, j+1, edgek+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<<e)) {
// is the vertex already calculated?
if(*eVert[ e ] < 0) {
// interpolate edge
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1
const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2
const float valp1 = value[ e1 ]; // scalar field val 1
const float valp2 = value[ e2 ]; // scalar field val 2
const float mu = (mIsoValue - valp1) / (valp2 - valp1);
// init isolevel vertex
ilv.v = p1 + (p2-p1)*mu;
ilv.n = getNormal( i+cubieOffsetX[e1], j+cubieOffsetY[e1], k+cubieOffsetZ[e1]) * (1.0-mu) +
getNormal( i+cubieOffsetX[e2], j+cubieOffsetY[e2], k+cubieOffsetZ[e2]) * ( mu) ;
mPoints.push_back( ilv );
triIndices[e] = (mPoints.size()-1);
// store vertex
*eVert[ e ] = triIndices[e];
} else {
// retrieve from vert array
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
if( (i<coAdd+mCutoff) || (j<coAdd+mCutoff) ||
((mCutoff>0) && (k<coAdd)) ||// bottom layer
(i>mSizex-2-coAdd-mCutoff) ||
(j>mSizey-2-coAdd-mCutoff) ) {
if(mCutArray) {
if(k < mCutArray[j*this->mSizex+i]) continue;
} else { continue; }
}
// Create the triangles...
for(int e=0; mcTriTable[cubeIndex][e]!=-1; e+=3) {
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<ParticleObject>::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<ParticleObject>::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:"<<pInUse<<", Subdivs:"<<mSubdivs, 9);
pz = mStart[2]-(double)(0.*gsz)-0.5*orgGsz;
for(int ok=1;ok<(mSizez-2)*mSubdivs;ok++) {
pz += gsz;
const int k = ok/mSubdivs;
if(k<=0) continue; // skip zero plane
for(int j=1;j<(mSizey-2);j++) {
for(int i=1;i<(mSizex-2);i++) {
orgval[0] = *getData(i ,j ,k );
orgval[1] = *getData(i+1,j ,k );
orgval[2] = *getData(i+1,j+1,k ); // with subdivs
orgval[3] = *getData(i ,j+1,k );
orgval[4] = *getData(i ,j ,k+1);
orgval[5] = *getData(i+1,j ,k+1);
orgval[6] = *getData(i+1,j+1,k+1); // with subdivs
orgval[7] = *getData(i ,j+1,k+1);
// prebuild subsampled array slice
const int sdkOffset = ok-k*mSubdivs;
for(int sdk=0; sdk<2; sdk++)
for(int sdj=0; sdj<mSubdivs+1; sdj++)
for(int sdi=0; sdi<mSubdivs+1; sdi++) {
subdAr[sdk][sdj][sdi] = ISOTRILININT(sdi*subdfac, sdj*subdfac, (sdkOffset+sdk)*subdfac);
}
const int poDistOffset=2;
for(int pok=-poDistOffset; pok<1+poDistOffset; pok++) {
if(k+pok<0) continue;
if(k+pok>=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<ParticleObject>::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(pfLen<minPfLen) pfLen = minPfLen;
//errMsg("ISOPPP"," at "<<PRINT_IJK<<" pp"<<ppos<<" sp"<<PRINT_VEC(spi,spj,spk)<<" pflen"<<pfLen );
//errMsg("ISOPPP"," subdfac="<<subdfac<<" size"<<p->getSize()<<" ps"<<mPartSize );
const int icellpsize = (int)(1.*pfLen*(gfxReal)mSubdivs)+1;
for(int swk=-icellpsize; swk<=icellpsize; swk++) {
if(spk+swk< 0) { continue; }
if(spk+swk> 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(len<pfLen) {
isoadd = baseIsoVal*1.;
} else {
// falloff linear with pfLen (kernel size=2pfLen
isoadd = baseIsoVal*(1. - (len-pfLen)/(pfLen));
}
if(isoadd<0.) { continue; }
//errMsg("ISOPPP"," at "<<PRINT_IJK<<" sp"<<PRINT_VEC(spi+swi,spj+swj,spk+swk)<<" cellp"<<cellp<<" pp"<<ppos << " l"<< len<< " add"<< isoadd);
const gfxReal arval = subdAr[spk+swk][spj+swj][spi+swi];
if(arval>1.) { 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<mSubdivs;sj++) {
py += gsy;
px = mStart[0]+(((double)i-0.5)*orgGsx)-gsx;
for(int si=0;si<mSubdivs;si++) {
px += gsx;
value[0] = subdAr[0+0][sj+0][si+0];
value[1] = subdAr[0+0][sj+0][si+1];
value[2] = subdAr[0+0][sj+1][si+1];
value[3] = subdAr[0+0][sj+1][si+0];
value[4] = subdAr[0+1][sj+0][si+0];
value[5] = subdAr[0+1][sj+0][si+1];
value[6] = subdAr[0+1][sj+1][si+1];
value[7] = subdAr[0+1][sj+1][si+0];
// check intersections of isosurface with edges, and calculate cubie index
cubeIndex = 0;
if (value[0] < mIsoValue) cubeIndex |= 1;
if (value[1] < mIsoValue) cubeIndex |= 2; // with subdivs
if (value[2] < mIsoValue) cubeIndex |= 4;
if (value[3] < mIsoValue) cubeIndex |= 8;
if (value[4] < mIsoValue) cubeIndex |= 16;
if (value[5] < mIsoValue) cubeIndex |= 32; // with subdivs
if (value[6] < mIsoValue) cubeIndex |= 64;
if (value[7] < mIsoValue) cubeIndex |= 128;
if (mcEdgeTable[cubeIndex] > 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<<e)) {
// is the vertex already calculated?
if(*eVert[ e ] < 0) {
// interpolate edge
const int e1 = mcEdges[e*2 ];
const int e2 = mcEdges[e*2+1];
const ntlVec3Gfx p1 = pos[ e1 ]; // scalar field pos 1
const ntlVec3Gfx p2 = pos[ e2 ]; // scalar field pos 2
const float valp1 = value[ e1 ]; // scalar field val 1
const float valp2 = value[ e2 ]; // scalar field val 2
const float mu = (mIsoValue - valp1) / (valp2 - valp1);
// init isolevel vertex
ilv.v = p1 + (p2-p1)*mu; // with subdivs
mPoints.push_back( ilv );
triIndices[e] = (mPoints.size()-1);
// store vertex
*eVert[ e ] = triIndices[e];
} else {
// retrieve from vert array
triIndices[e] = *eVert[ e ];
}
} // along all edges
}
// removed cutoff treatment...
// 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] ] ); // with subdivs
mIndices.push_back( triIndices[ mcTriTable[cubeIndex][e+2] ] );
//errMsg("TTT"," i1"<<mIndices[mIndices.size()-3]<<" "<< " i2"<<mIndices[mIndices.size()-2]<<" "<< " i3"<<mIndices[mIndices.size()-1]<<" "<< mIndices.size() );
}
} // triangles in edge table?
}//si
}// sj
}//i
}// j
// copy edge arrays
for(int j=0;j<(mSizey-0)*mSubdivs;j++)
for(int i=0;i<(mSizex-0)*mSubdivs;i++) {
//int edgek = 0;
const int dst = EDGEAR_INDEX( 0, 0, 0, i,j);
const int src = EDGEAR_INDEX( 0, 0, 1, i,j);
mpEdgeVerticesX[ dst ] = mpEdgeVerticesX[ src ];
mpEdgeVerticesY[ dst ] = mpEdgeVerticesY[ src ]; // with subdivs
mpEdgeVerticesZ[ dst ] = mpEdgeVerticesZ[ src ];
mpEdgeVerticesX[ src ]=-1;
mpEdgeVerticesY[ src ]=-1; // with subdivs
mpEdgeVerticesZ[ src ]=-1;
}
// */
} // ok, k subdiv loop
//delete [] subdAr;
delete [] arppnt;
computeNormals();
} // with subdivs
// perform smoothing
float smoSubdfac = 1.;
if(mSubdivs>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("<<mSmoothSurface<<","<<mSmoothNormals<<")" , 10 );
debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< getTimeString(tritimeend-tritimestart)<<", S("<<mSmoothSurface<<","<<mSmoothNormals<<"),"<<
" verts:"<<mPoints.size()<<" tris:"<<(mIndices.size()/3)<<" subdivs:"<<mSubdivs
, 10 );
if(mpIsoParts) debMsgStd("IsoSurface::triangulate",DM_MSG,"parts:"<<mpIsoParts->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"<<srcp<<" c"<<mSCcenter<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" ");
//invsigma2 *= (rd*rd-mSCrad1);
//flt = ntlVec3Gfx(100); return 1;
return 0;
} else if(rd > 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"<<srcp<<" rd:"<<rd<<" r1:"<<mSCrad1<<" r2:"<<mSCrad2<<" org:"<<org<<" is:"<<invsigma2);
//invsigma2 *= (rd*rd-mSCrad1);
//return 0;
} else {
}
}
flt = ntlVec3Gfx(0.0);
flt += *(field+pointerScale*src) *pointareas[src];
float sum_w = pointareas[src];
//const ntlVec3Gfx &nv = mPoints[src].n;
target = ntlVec3Gfx(0.0);
target += *(field+pointerScale*src) *pointareas[src];
float smstrSum = pointareas[src];
int flag = mFlagCnt;
mFlagCnt++;
flags[src] = flag;
//vector<int> 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:"<<sigma); // DEBUG
}
void IsoSurface::smoothNormals(float sigma)
{
// only smoothen the normals
void IsoSurface::smoothNormals(float sigma) {
// reuse from smoothSurface
if(neighbors.size() != mPoints.size()) {
// need neighbor
@@ -781,10 +1123,8 @@ void IsoSurface::smoothNormals(float sigma)
} else { nflt[i]=mPoints[i].n; }
}
// copy back
// copy results
for (int i = 0; i < nv; i++) { mPoints[i].n = nflt[i]; }
//errMsg("SMNRMLS","done v:"<<sigma); // DEBUG
}

27
intern/elbeem/intern/isosurface.h
View File

@@ -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
*
* Marching Cubes "displayer"
*
@@ -18,6 +18,8 @@
/* access some 3d array */
#define ISOLEVEL_INDEX(ii,ij,ik) ((mSizex*mSizey*(ik))+(mSizex*(ij))+((ii)))
class ParticleTracer;
/* struct for a small cube in the scalar field */
typedef struct {
ntlVec3Gfx pos[8];
@@ -53,6 +55,19 @@ class IsoSurface :
/*! triangulate the scalar field given by pointer*/
void triangulate( void );
/*! set particle pointer */
void setParticles(ParticleTracer *pnt,float psize){ mpIsoParts = pnt; mPartSize=psize; };
/*! set # of subdivisions, this has to be done before init! */
void setSubdivs(int s) {
if(mInitDone) errFatal("IsoSurface::setSubdivs","Changing subdivs after init!", SIMWORLD_INITERROR);
if(s<1) s=1; if(s>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<IsoLevelVertex> 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<int> 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:

4
intern/elbeem/intern/loop_tools.h
View File

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

69
intern/elbeem/intern/ntl_blenderdumper.cpp
View File

@@ -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 <fstream>
#include <sys/types.h>
@@ -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<ntlVec3Gfx> Vertices;
vector<ntlVec3Gfx> VertNormals;
/* init geometry array, first all standard objects */
// check geo objects
int idCnt = 0; // give IDs to objects
for (vector<ntlGeometryClass*>::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 "<<geoobj->getName(), 8 );
// normal geom. objects -> ignore
}
if(tid & GEOCLASSTID_SHADER) {
ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter); //dynamic_cast<ntlGeometryShader*>(*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 '"<<boutfilename<<"' ");
return 1; }
//! current transform matrix
// dont transform velocity output, this is handled in blender
// current transform matrix
ntlMatrix4x4<gfxReal> *trafo;
trafo = lbm->getDomainTrafo();
if(trafo) {
//trafo->initArrayCheck(ettings->surfaceTrafo);
//errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) );
// transform into source space
for(size_t i=0; i<Vertices.size(); i++) {
Vertices[i] = (*trafo) * Vertices[i];
}
}
// rotate vertnormals
ntlMatrix4x4<gfxReal> 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; i<Vertices.size(); i++) {
VertNormals[i] = rottrafo * VertNormals[i];
normalize(VertNormals[i]); // remove scaling etc.
}
}
// write to file
int numVerts;
if(sizeof(numVerts)!=4) { errMsg("ntlBlenderDumper::renderScene","Invalid int size"); return 1; }
@@ -215,18 +236,20 @@ int ntlBlenderDumper::renderScene( void )
if(numGMs>0) {
if(debugOut) debMsgStd("ntlBlenderDumper::renderScene",DM_MSG,"Objects dumped: "<<numGMs, 10);
} else {
errFatal("ntlBlenderDumper::renderScene","No objects to dump! Aborting...",SIMWORLD_INITERROR);
return 1;
if((glob_mpactive) && (glob_mpindex>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 #"<<nrStr<<" dump time: "<< getTimeString(stopTime-startTime) <<" ", 10);
// still render for preview...
debugOut = false; // DEBUG!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(debugOut) {
if(debugRender) {
debMsgStd("ntlBlenderDumper::renderScene",DM_NOTIFY,"Performing preliminary render", 1);
ntlWorld::renderScene(); }
else {

2
intern/elbeem/intern/ntl_blenderdumper.h
View File

@@ -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
*
* Replaces std. raytracer, and only dumps time dep. objects to disc, header
*

2
intern/elbeem/intern/ntl_bsptree.cpp
View File

@@ -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
*
* Tree container for fast triangle intersects
*

2
intern/elbeem/intern/ntl_bsptree.h
View File

@@ -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
*
* Tree container for fast triangle intersects
*

9
intern/elbeem/intern/ntl_geometryclass.h
View File

@@ -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
*
* Base class for geometry shaders and objects
*
@@ -31,6 +31,7 @@ class ntlGeometryClass
mObjectId(-1), mpAttrs( NULL ), mGeoInitId(-1)
{
mpAttrs = new AttributeList("objAttrs");
mpSwsAttrs = new AttributeList("swsAttrs");
};
//! Default destructor
@@ -58,6 +59,10 @@ class ntlGeometryClass
/*! Returns the attribute list pointer */
inline AttributeList *getAttributeList() { return mpAttrs; }
/*! Get/Sets the attribute list pointer */
inline void setSwsAttributeList(AttributeList *set) { mpSwsAttrs=set; }
inline AttributeList *getSwsAttributeList() { return mpSwsAttrs; }
/*! for easy GUI detection get start of axis aligned bounding box, return NULL of no BB */
virtual inline ntlVec3Gfx *getBBStart() { return NULL; }
virtual inline ntlVec3Gfx *getBBEnd() { return NULL; }
@@ -93,6 +98,8 @@ class ntlGeometryClass
/*! configuration attributes */
AttributeList *mpAttrs;
/*! sws configuration attributes */
AttributeList *mpSwsAttrs;
/* fluid init data */
/*! id of fluid init (is used in solver initialization), additional data stored only for objects */

2
intern/elbeem/intern/ntl_geometrymodel.cpp
View File

@@ -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
*
* A simple box object
*

2
intern/elbeem/intern/ntl_geometrymodel.h
View File

@@ -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
*
* A model laoded from Wavefront .obj file
*

51
intern/elbeem/intern/ntl_geometryobject.cpp
View File

@@ -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
*
* a geometry object
* all other geometry objects are derived from this one
@@ -16,6 +16,8 @@
// for FGI
#include "elbeem.h"
#define TRI_UVOFFSET (1./4.)
//#define TRI_UVOFFSET (1./3.)
/*****************************************************************************/
@@ -29,6 +31,7 @@ ntlGeometryObject::ntlGeometryObject() :
mInitialVelocity(0.0), mcInitialVelocity(0.0), mLocalCoordInivel(false),
mGeoInitIntersect(false),
mGeoPartSlipValue(0.0),
mcGeoImpactFactor(1.),
mVolumeInit(VOLUMEINIT_VOLUME),
mInitialPos(0.),
mcTrans(0.), mcRot(0.), mcScale(1.),
@@ -62,6 +65,7 @@ bool ntlGeometryObject::checkIsAnimated() {
|| (mcRot.accessValues().size()>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="<<getName()<<" debug: trans:"<<mcTrans.accessValues().size()<<" rot:"<<mcRot.accessValues().size()<<" scale:"<<mcScale.accessValues().size()<<" active:"<<mcGeoActive.accessValues().size()<<" inivel:"<<mcInitialVelocity.accessValues().size()<<". isani?"<<mIsAnimated ); // DEBUG
return mIsAnimated;
}
@@ -139,6 +144,13 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
mVolumeInit = mpAttrs->readInt("geoinit_volumeinit", mVolumeInit,"ntlGeometryObject", "mVolumeInit", false);
if((mVolumeInit<VOLUMEINIT_VOLUME)||(mVolumeInit>VOLUMEINIT_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<double>(geoactive); }
if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel<float>(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:"<<nTrans<<" nr:"<<nRot<<" ns:"<<nScale, 10);
debMsgStd("ntlGeometryObject::initChannels",DM_MSG,"na:"<<nAct<<" niv:"<<nIvel<<" ", 10); }
vector<ntlVec3Gfx> vals;
vector<double> valsd;
vector<float> valsfloat;
vector<double> 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<ntlVec3Gfx> *verts,
ntlMat4Gfx rotMat;
rotMat.initRotationXYZ(rot[0],rot[1],rot[2]);
pos += mInitialPos;
//errMsg("ntlGeometryObject::applyTransformation","obj="<<getName()<<" t"<<pos<<" r"<<rot<<" s"<<scale);
errMsg("ntlGeometryObject::applyTransformation","obj="<<getName()<<" t"<<pos<<" r"<<rot<<" s"<<scale);
for(int i=vstart; i<vend; i++) {
(*verts)[i] *= scale;
(*verts)[i] = rotMat * (*verts)[i];
@@ -397,7 +409,7 @@ void ntlGeometryObject::applyTransformation(double t, vector<ntlVec3Gfx> *verts,
}
} else {
// not animated, cached points were already returned
//errMsg ("ntlGeometryObject::applyTransformation","Object "<<getName()<<" used cached points ");
errMsg ("ntlGeometryObject::applyTransformation","Object "<<getName()<<" used cached points ");
}
}
@@ -473,8 +485,8 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
if(divs1+divs2 > 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"<<i<<" p"<<p<<" n"<<n<<" ");
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","std"<<i<<" p"<<p<<" n"<<n<<" ");
}
// init points & refine...
for(size_t i=0; i<triangles.size(); i++) {
@@ -515,12 +527,12 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
if(discardInflowBack) {
if(dot(mInitialVelocity,trinorm)<0.0) continue;
}
//errMsg("ntlGeometryObject::initMovingPoints","Tri1 "<<vertices[trips[0]]<<","<<vertices[trips[1]]<<","<<vertices[trips[2]]<<" "<<divs1<<","<<divs2 );
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Tri1 "<<vertices[trips[0]]<<","<<vertices[trips[1]]<<","<<vertices[trips[2]]<<" "<<divs1<<","<<divs2 );
if(divs1+divs2 > 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 "<<srcvertices[srctrips[0]]<<","<<srcvertices[srctrips[1]]<<","<<srcvertices[srctrips[2]]<<" "<<divs1<<","<<divs2 );
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 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="<<t1<<" t2="<<t2<<" p1="<<verts1[0]<<" p2="<<verts2[0]<<" v="<<vel);
//errMsg("ntlGeometryObject::calculateMaxVel","t1="<<t1<<" t2="<<t2<<" p1="<<verts1[0]<<" p2="<<verts2[0]<<" v="<<vel);
return vel;
}
@@ -758,7 +770,6 @@ ntlVec3Gfx ntlGeometryObject::getTranslation(double t) {
}
/*! get active flag time t*/
float ntlGeometryObject::getGeoActive(double t) {
//float act = mcGeoActive.getConstant(t);
float act = mcGeoActive.get(t); // if <= 0.0 -> off
return act;
}

10
intern/elbeem/intern/ntl_geometryobject.h
View File

@@ -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
*
* a geometry object
* all other geometry objects are derived from this one
@@ -82,6 +82,10 @@ class ntlGeometryObject : public ntlGeometryClass
inline float getGeoPartSlipValue() const { return mGeoPartSlipValue; }
inline void setGeoPartSlipValue(float set) { mGeoPartSlipValue=set; }
/*! Set/get the impact corr factor channel */
inline float getGeoImpactFactor(double t) { return mcGeoImpactFactor.get(t); }
inline void setGeoImpactFactor(float set) { mcGeoImpactFactor = AnimChannel<float>(set); }
/*! Set/get the part slip value*/
inline int getVolumeInit() const { return mVolumeInit; }
inline void setVolumeInit(int set) { mVolumeInit=set; }
@@ -170,6 +174,8 @@ class ntlGeometryObject : public ntlGeometryClass
bool mGeoInitIntersect;
/*! part slip bc value */
float mGeoPartSlipValue;
/*! obstacle impact correction factor */
AnimChannel<float> mcGeoImpactFactor;
/*! only init as thin object, dont fill? */
int mVolumeInit;
@@ -195,7 +201,7 @@ class ntlGeometryObject : public ntlGeometryClass
int mMaxMovPnt;
/*! animated channels for in/outflow on/off */
AnimChannel<double> mcGeoActive;
AnimChannel<float> mcGeoActive;
public:

2
intern/elbeem/intern/ntl_geometryshader.h
View File

@@ -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
*
* Interface for a geometry shader
*

2
intern/elbeem/intern/ntl_lighting.cpp
View File

@@ -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
*
* a light object
*

2
intern/elbeem/intern/ntl_lighting.h
View File

@@ -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
*
* a light object
* default omni light implementation

81
intern/elbeem/intern/ntl_matrices.h
View File

@@ -2,7 +2,7 @@
/******************************************************************************
*
* El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
* Copyright 2003,2004 Nils Thuerey
* Copyright 2003-2006 Nils Thuerey
*
* Basic matrix utility include file
*
@@ -84,6 +84,9 @@ public:
//! from 16 value array (init id if all 0)
inline void initArrayCheck(Scalar *array);
//! decompose matrix
void decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear);
//! public to avoid [][] operators
Scalar value[4][4]; //< Storage of vector values
@@ -695,6 +698,82 @@ ntlMatrix4x4<Scalar>::initArrayCheck(Scalar *array)
if(allZero) this->initId();
}
//! decompose matrix
template<class Scalar>
void
ntlMatrix4x4<Scalar>::decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear) {
ntlVec3Gfx row[3],temp;
for(int i = 0; i < 3; i++) {
trans[i] = this->value[3][i];
}
for(int i = 0; i < 3; i++) {
row[i][0] = this->value[i][0];
row[i][1] = this->value[i][1];
row[i][2] = this->value[i][2];
}
scale[0] = norm(row[0]);
normalize (row[0]);
shear[0] = dot(row[0], row[1]);
row[1][0] = row[1][0] - shear[0]*row[0][0];
row[1][1] = row[1][1] - shear[0]*row[0][1];
row[1][2] = row[1][2] - shear[0]*row[0][2];
scale[1] = norm(row[1]);
normalize (row[1]);
if(scale[1] != 0.0)
shear[0] /= scale[1];
shear[1] = dot(row[0], row[2]);
row[2][0] = row[2][0] - shear[1]*row[0][0];
row[2][1] = row[2][1] - shear[1]*row[0][1];
row[2][2] = row[2][2] - shear[1]*row[0][2];
shear[2] = dot(row[1], row[2]);
row[2][0] = row[2][0] - shear[2]*row[1][0];
row[2][1] = row[2][1] - shear[2]*row[1][1];
row[2][2] = row[2][2] - shear[2]*row[1][2];
scale[2] = norm(row[2]);
normalize (row[2]);
if(scale[2] != 0.0) {
shear[1] /= scale[2];
shear[2] /= scale[2];
}
temp = cross(row[1], row[2]);
if(dot(row[0], temp) < 0.0) {
for(int i = 0; i < 3; i++) {
scale[i] *= -1.0;
row[i][0] *= -1.0;
row[i][1] *= -1.0;
row[i][2] *= -1.0;
}
}
if(row[0][2] < -1.0) row[0][2] = -1.0;
if(row[0][2] > +1.0) row[0][2] = +1.0;
rot[1] = asin(-row[0][2]);
if(fabs(cos(rot[1])) > VECTOR_EPSILON) {
rot[0] = atan2 (row[1][2], row[2][2]);
rot[2] = atan2 (row[0][1], row[0][0]);
}
else {
rot[0] = atan2 (row[1][0], row[1][1]);
rot[2] = 0.0;
}
rot[0] = (180.0/M_PI)*rot[0];
rot[1] = (180.0/M_PI)*rot[1];
rot[2] = (180.0/M_PI)*rot[2];
}
#define NTL_MATRICES_H
#endif

2
intern/elbeem/intern/ntl_ray.cpp
View File

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

2
intern/elbeem/intern/ntl_ray.h
View File

@@ -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
*
* ray class
*

2
intern/elbeem/intern/ntl_vector3dim.h
View File

@@ -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 vector class used everywhere, either blitz or inlined GRAPA class
*

26
intern/elbeem/intern/ntl_world.cpp
View File

@@ -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
*
@@ -421,6 +421,7 @@ int ntlWorld::advanceSims(int framenum)
//gstate = getGlobalBakeState();
//if(gstate<0) { allPanic = true; done = true; } // this means abort... cause panic
//#endif // ELBEEM_BLENDER==1
myTime_t advsstart = getTime();
// step all the sims, and check for panic
debMsgStd("ntlWorld::advanceSims",DM_MSG, " sims "<<mpSims->size()<<" t"<<targetTime<<" done:"<<done<<" panic:"<<allPanic<<" gstate:"<<gstate, 10); // debug // timedebug
@@ -453,6 +454,9 @@ int ntlWorld::advanceSims(int framenum)
return 1;
}
myTime_t advsend = getTime();
debMsgStd("ntlWorld::advanceSims",DM_MSG,"Overall steps so far took:"<< getTimeString(advsend-advsstart)<<" for sim time "<<targetTime, 4);
// finish step
for(size_t i=0;i<mpSims->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() <<"_"<<glob_mpindex<<"_" << nrStr << ".png"; /// DEBUG!
} else {
// ORG
outfn_conv << glob->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(),

2
intern/elbeem/intern/ntl_world.h
View File

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

16
intern/elbeem/intern/parametrizer.cpp
View File

@@ -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="<<mCellSize ,10);
@@ -424,8 +426,8 @@ errMsg("Warning","Used z-dir for gstar!");
if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," gravity force = "<<PRINT_NTLVEC(mcGravity.get(time))<<", scaled with "<<forceFactor<<" to "<<calculateGravity(time),1);
}
mExtent = ParamVec( mCellSize*mSizex, mCellSize*mSizey, mCellSize*mSizez );
if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," domain extent = "<<PRINT_NTLVEC(mExtent)<<"m , gs:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<" cs:"<<mCellSize,1);
//mExtent = ParamVec( mCellSize*mSizex, mCellSize*mSizey, mCellSize*mSizez );
//if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," domain extent = "<<PRINT_NTLVEC(mExtent)<<"m , gs:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<" cs:"<<mCellSize,1);
if(!checkSeenValues(PARAM_ANIFRAMETIME)) {
errFatal("Parametrizer::calculateAllMissingValues"," Warning no ani frame time given!", SIMWORLD_INITERROR);

6
intern/elbeem/intern/parametrizer.h
View File

@@ -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
*
@@ -104,7 +104,7 @@ class Parametrizer {
/*! get start time of animation */
int calculateAniStart( void );
/*! get extent of the domain = (1,1,1) if parametrizer not used, (x,y,z) [m] otherwise */
ParamVec calculateExtent( void );
//ParamVec calculateExtent( void );
/*! get (scaled) surface tension */
ParamFloat calculateSurfaceTension( void );
/*! get time step size for lbm (equals mTimeFactor) */
@@ -271,7 +271,7 @@ class Parametrizer {
ParamFloat mAniStart;
/*! extent of the domain in meters */
ParamVec mExtent;
//ParamVec mExtent;
/*! fluid density [kg/m^3], default 1.0 g/cm^3 */
ParamFloat mDensity;

22
intern/elbeem/intern/particletracer.cpp
View File

@@ -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
*
@@ -128,12 +128,14 @@ void ParticleTracer::adaptPartTimestep(float factor) {
/******************************************************************************
* add a particle at this position
*****************************************************************************/
void ParticleTracer::addParticle(float x, float y, float z)
{
void ParticleTracer::addParticle(float x, float y, float z) {
ntlVec3Gfx p(x,y,z);
ParticleObject part( p );
mParts.push_back( part );
}
void ParticleTracer::addFullParticle(ParticleObject &np) {
mParts.push_back( np );
}
void ParticleTracer::cleanup() {
@@ -150,6 +152,9 @@ void ParticleTracer::cleanup() {
}
}
extern bool glob_mpactive;
extern int glob_mpindex,glob_mpnum;
/******************************************************************************
*! dump particles if desired
*****************************************************************************/
@@ -161,8 +166,13 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
(mDumpParts>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"<<glob_mpindex; }
}
boutfilename << ".gz";
debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
//debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: partgeodeb sim:"<<mSimStart<<","<<mSimEnd<<" geosize:"<<mStart<<","<<mEnd,2 );
// output to zipped file
gzFile gzf;
@@ -211,7 +221,9 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
void ParticleTracer::checkDumpTextPositions(double simtime) {
// dfor partial & full dump
errMsg("ParticleTracer::checkDumpTextPositions","t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval);
if(mDumpTextInterval>0.) {
debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval,7);
}
if((mDumpTextInterval>0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) {
// dump to binary file

21
intern/elbeem/intern/particletracer.h
View File

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

115
intern/elbeem/intern/simulation_object.cpp
View File

@@ -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:"<<mpLbm->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: "<<dinitType<<" ps="<<mpElbeemSettings->obstaclePartslip<<" vv"<<mpElbeemSettings->generateVertexVectors<<","<<mpLbm->getDumpVelocities(), 9 );
debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<<dinitType<<" ps="<<mpElbeemSettings->domainobsPartslip<<" vv"<<mpElbeemSettings->generateVertexVectors<<","<<mpLbm->getDumpVelocities(), 9 );
debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<<mpLbm->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; j<jmax ; j++) flagCount[j] = 0;
int diffInits = 0;
LbmSolverInterface::CellIdentifier cid = mpLbm->getFirstCell();
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; j<jmax ; j++) flagCount[j] = 0;
int diffInits = 0;
LbmSolverInterface::CellIdentifier cid = mpLbm->getFirstCell();
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; j<jmax ; j++) {
if( flag&(1<<j) ) flagCount[j]++;
}
totalCells++;
}
for(int j=0; j<jmax ; j++) {
if( flag&(1<<j) ) flagCount[j]++;
}
totalCells++;
}
mpLbm->deleteCellIterator( &cid );
mpLbm->deleteCellIterator( &cid );
char charNl = '\n';
debugOutNnl("SimulationObject::initializeLbmSimulation celltype stats: " <<charNl, 5);
debugOutNnl("no. of cells = "<<totalCells<<", "<<charNl ,5);
for(int j=0; j<jmax ; j++) {
std::ostringstream out;
if(flagCount[j]>0) {
out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1<<j) ) <<", " << charNl;
char charNl = '\n';
debugOutNnl("SimulationObject::initializeLbmSimulation celltype stats: " <<charNl, 5);
debugOutNnl("no. of cells = "<<totalCells<<", "<<charNl ,5);
for(int j=0; j<jmax ; j++) {
std::ostringstream out;
if(flagCount[j]>0) {
out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1<<j) ) <<", " << charNl;
debugOutNnl(out.str(), 5);
}
}
// compute dist. of empty/bnd - fluid - if
// cfEmpty = (1<<0), cfBnd = (1<< 2), cfFluid = (1<<10), cfInter = (1<<11),
if(1){
std::ostringstream out;
out.precision(2); out.width(4);
int totNum = flagCount[1]+flagCount[2]+flagCount[7]+flagCount[8];
double ebFrac = (double)(flagCount[1]+flagCount[2]) / totNum;
double flFrac = (double)(flagCount[7]) / totNum;
double ifFrac = (double)(flagCount[8]) / totNum;
//???
out<<"\tFractions: [empty/bnd - fluid - interface - ext. if] = [" << ebFrac<<" - " << flFrac<<" - " << ifFrac<<"] "<< charNl;
if(diffInits > 0) {
debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="<<diffInits<<"!" , 5);
}
debugOutNnl(out.str(), 5);
}
}
// compute dist. of empty/bnd - fluid - if
// cfEmpty = (1<<0), cfBnd = (1<< 2), cfFluid = (1<<10), cfInter = (1<<11),
if(1){
std::ostringstream out;
out.precision(2); out.width(4);
int totNum = flagCount[1]+flagCount[2]+flagCount[7]+flagCount[8];
double ebFrac = (double)(flagCount[1]+flagCount[2]) / totNum;
double flFrac = (double)(flagCount[7]) / totNum;
double ifFrac = (double)(flagCount[8]) / totNum;
//???
out<<"\tFractions: [empty/bnd - fluid - interface - ext. if] = [" << ebFrac<<" - " << flFrac<<" - " << ifFrac<<"] "<< charNl;
if(diffInits > 0) {
debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="<<diffInits<<"!" , 5);
}
debugOutNnl(out.str(), 5);
}
} // cellstats
// might be modified by mpLbm
mpParts->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...

2
intern/elbeem/intern/simulation_object.h
View File

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

146
intern/elbeem/intern/solver_adap.cpp
View File

@@ -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:"<<glob_mpactive<<","<<glob_mpindex<<"/"<<glob_mpnum<<" step:"<<mLevel[mMaxRefine].lsteps<<" skipping tadap...",8);
return;
}
nextmax = mrInitTadap(nextmax);
}
#endif // LBM_INCLUDE_TESTSOLVERS==1
LbmFloat newdt = mpParam->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(nextmax<allowMax*reduceFac) {
newdt = this->mpParam->getTimestep() / reduceFac;
newdt = mpParam->getTimestep() / reduceFac;
}
} // newtr
//errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< this->mpParam->getSimulationMaxSpeed() );
//errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< mpParam->getSimulationMaxSpeed() );
bool minCutoff = false;
LbmFloat desireddt = newdt;
if(newdt>this->mpParam->getMaxTimestep()){ newdt = this->mpParam->getMaxTimestep(); }
if(newdt<this->mpParam->getMinTimestep()){
newdt = this->mpParam->getMinTimestep();
if(newdt>mpParam->getMaxTimestep()){ newdt = mpParam->getMaxTimestep(); }
if(newdt<mpParam->getMinTimestep()){
newdt = mpParam->getMinTimestep();
if(nextmax>allowMax/reduceFac){ minCutoff=true; } // only if really large vels...
}
LbmFloat dtdiff = fabs(newdt - this->mpParam->getTimestep());
if(!this->mSilent) {
LbmFloat dtdiff = fabs(newdt - mpParam->getTimestep());
if(!mSilent) {
debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt
<<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff
<<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
<<" max"<<mpParam->getMaxTimestep()<<" min"<<mpParam->getMinTimestep()<<" diff"<<dtdiff
<<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep)<<
" olddt"<<levOldStepsize[mMaxRefine]<<" redlock"<<mTimestepReduceLock
, 10); }
// in range, and more than X% change?
//if( newdt < this->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... "<<mTimestepReduceLock<<" ",10);
debMsgDirect("D");
} else {
this->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="<<newdt<<" old="<<this->mpParam->getTimestep()
<<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<<
"rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: new="<<newdt<<" old="<<mpParam->getTimestep()
<<" maxSpeed:"<<mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<mStepCnt, 10 );
//debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: "<< "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
}
} // really change dt
}
@@ -1009,17 +1023,17 @@ void LbmFsgrSolver::adaptTimestep() {
/*const int tadtogInter = 100;
const double tadtogSwitch = 0.66;
errMsg("TIMESWITCHTOGGLETEST","warning enabled "<< tadtogSwitch<<","<<tadtogSwitch<<" !!!!!!!!!!!!!!!!!!!");
if( ((this->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()<mMinTimestep) mMinTimestep = this->mpParam->getTimestep();
if(this->mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = this->mpParam->getTimestep();
if(mpParam->getTimestep()<mMinTimestep) mMinTimestep = mpParam->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: "<<lev<<" Timestep change: "<<
if(!mSilent) {
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep chlevel: "<<
" scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10);
}
if(lev!=mMaxRefine) coarseCalculateFluxareas(lev);
@@ -1079,7 +1095,7 @@ void LbmFsgrSolver::adaptTimestep() {
// init empty zero vel interface cell...
initVelocityCell(lev, i,j,k, CFInter, 1.0, 0.01, LbmVec(0.) );
} else {// */
for(int l=0; l<this->cDfNum; l++) {
for(int l=0; l<cDfNum; l++) {
QCELL(lev, i, j, k, workSet, l) = QCELL(lev, i, j, k, workSet, l)* scaleFac;
}
//} // ok
@@ -1102,12 +1118,12 @@ void LbmFsgrSolver::adaptTimestep() {
LbmVec velOld;
LbmFloat rho, ux,uy,uz;
rho=0.0; ux = uy = uz = 0.0;
for(int l=0; l<this->cDfNum; l++) {
for(int l=0; l<cDfNum; l++) {
LbmFloat m = QCELL(lev, i, j, k, workSet, l);
rho += m;
ux += (this->dfDvecX[l]*m);
uy += (this->dfDvecY[l]*m);
uz += (this->dfDvecZ[l]*m);
ux += (dfDvecX[l]*m);
uy += (dfDvecY[l]*m);
uz += (dfDvecZ[l]*m);
}
rhoOld = rho;
velOld = LbmVec(ux,uy,uz);
@@ -1118,21 +1134,21 @@ void LbmFsgrSolver::adaptTimestep() {
LbmFloat df[LBM_DFNUM];
LbmFloat feqOld[LBM_DFNUM];
LbmFloat feqNew[LBM_DFNUM];
for(int l=0; l<this->cDfNum; l++) {
feqOld[l] = this->getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
feqNew[l] = this->getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
for(int l=0; l<cDfNum; l++) {
feqOld[l] = getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
feqNew[l] = getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
df[l] = QCELL(lev, i,j,k,workSet, l);
}
const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feqOld);
const LbmFloat oldOmega = this->getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
const LbmFloat newOmega = this->getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
const LbmFloat Qo = getLesNoneqTensorCoeff(df,feqOld);
const LbmFloat oldOmega = getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
const LbmFloat newOmega = getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
//newOmega = mLevel[lev].omega; // FIXME debug test
//LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
const LbmFloat dfScale = (newSteptime/newOmega)/(levOldStepsize[lev]/oldOmega);
//dfScale = dfScaleFac/newOmega;
for(int l=0; l<this->cDfNum; l++) {
for(int l=0; l<cDfNum; l++) {
// org scaling
//df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important
// new scaling
@@ -1176,22 +1192,22 @@ void LbmFsgrSolver::adaptTimestep() {
} // lev
if(!this->mSilent) {
debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<this->mStepCnt<<
if(!mSilent) {
debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<mStepCnt<<
" no"<<mTimeSwitchCounts<<" maxdt"<<mMaxTimestep<<
" mindt"<<mMinTimestep<<" currdt"<<mLevel[mMaxRefine].timestep, 10);
debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE masst:"<<massTNew<<","<<massTOld<<" org:"<<mCurrentMass<<"; "<<
" volt:"<<volTNew<<","<<volTOld<<" org:"<<mCurrentVolume, 10);
} else {
debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep change by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
<<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<this->mOmega<<" gStar:"<<this->mpParam->getCurrentGStar() , 10);
debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep changed by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
<<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<mOmega<<" gStar:"<<mpParam->getCurrentGStar()<<", step:"<<mStepCnt , 10);
}
} // rescale?
//NEWDEBCHECK("tt2");
//errMsg("adaptTimestep","Warning - brute force rescale off!"); minCutoff = false; // DEBUG
if(minCutoff) {
errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<this->mName<<" step:"<<this->mStepCnt<<" newdt="<<desireddt<<" mindt="<<this->mpParam->getMinTimestep()<<") " );
errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<mName<<" step:"<<mStepCnt<<" newdt="<<desireddt<<" mindt="<<mpParam->getMinTimestep()<<") " );
//brute force resacle all the time?
for(int lev=mMaxRefine; lev>=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; l<this->cDfNum; l++) {
for(int l=0; l<cDfNum; l++) {
LbmFloat m = QCELL(lev, i, j, k, workSet, l);
rho += m;
ux += (this->dfDvecX[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; l<this->cDfNum; l++) {
QCELL(lev, i, j, k, workSet, l) = this->getCollideEq(l, rho, ux,uy,uz); }
for(int l=0; l<cDfNum; l++) {
QCELL(lev, i, j, k, workSet, l) = getCollideEq(l, rho, ux,uy,uz); }
rescs++;
debMsgDirect("B");
}

78
intern/elbeem/intern/solver_class.h
View File

@@ -3,7 +3,7 @@
* El'Beem - the visual lattice boltzmann freesurface simulator
* 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
*
* Combined 2D/3D Lattice Boltzmann standard solver classes
*
@@ -259,13 +259,17 @@ class LbmFsgrSolver :
LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
LBM_INLINED void forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
LBM_INLINED void initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass);
//! interpolate velocity and density at a given position
void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz);
/*! perform a single LBM step */
void stepMain();
//! advance fine grid
void fineAdvance();
//! advance coarse grid
void coarseAdvance(int lev);
//! update flux area values on coarse grids
void coarseCalculateFluxareas(int lev);
// adaptively coarsen grid
bool adaptGrid(int lev);
@@ -278,6 +282,11 @@ class LbmFsgrSolver :
virtual int initParticles();
/*! move all particles */
virtual void advanceParticles();
/*! move a particle at a boundary */
void handleObstacleParticle(ParticleObject *p);
/*! check whether to add particle
bool checkAddParticle();
void performAddParticle();*/
/*! debug object display (used e.g. for preview surface) */
@@ -294,9 +303,6 @@ class LbmFsgrSolver :
//! for raytracing, preprocess
void prepareVisualization( void );
/*! type for cells */
//typedef typename this->LbmCell 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

555
intern/elbeem/intern/solver_init.cpp
View File

File diff suppressed because it is too large Load Diff

111
intern/elbeem/intern/solver_interface.cpp
View File

@@ -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
*
* Combined 2D/3D Lattice Boltzmann Interface Class
* contains stuff to be statically compiled
@@ -35,7 +35,7 @@ LbmSolverInterface::LbmSolverInterface() :
mBoundarySouth( (CellFlagType)(CFBnd) ),mBoundaryTop( (CellFlagType)(CFBnd) ),mBoundaryBottom( (CellFlagType)(CFBnd) ),
mInitDone( false ),
mInitDensityGradient( false ),
mpAttrs( NULL ), mpParam( NULL ), mpParticles(NULL),
mpSifAttrs( NULL ), mpSifSwsAttrs(NULL), mpParam( NULL ), mpParticles(NULL),
mNumParticlesLost(0),
mNumInvalidDfs(0), mNumFilledCells(0), mNumEmptiedCells(0), mNumUsedCells(0), mMLSUPS(0),
mDebugVelScale( 0.01 ), mNodeInfoString("+"),
@@ -50,17 +50,34 @@ LbmSolverInterface::LbmSolverInterface() :
mRefinementDesired(0),
mOutputSurfacePreview(0), mPreviewFactor(0.25),
mSmoothSurface(1.0), mSmoothNormals(1.0),
mPartGenProb(0.),
mIsoSubdivs(1), mPartGenProb(0.),
mDumpVelocities(false),
mMarkedCells(), mMarkedCellIndex(0),
mDomainBound("noslip"), mDomainPartSlipValue(0.1),
mTForceStrength(0.0), mFarFieldSize(0.), mCppfStage(0)
mFarFieldSize(0.),
mPartDropMassSub(0.1), // good default
mPartUsePhysModel(false),
mTForceStrength(0.0),
mCppfStage(0),
mDumpRawText(false),
mDumpRawBinary(false),
mDumpRawBinaryZip(true)
{
#if ELBEEM_PLUGIN==1
if(gDebugLevel<=1) mSilent = true;
if(gDebugLevel<=1) setSilent(true);
#endif
mpSimTrafo = new ntlMat4Gfx(0.0);
mpSimTrafo->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:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" " ,10);
const int debugGridsizeInit = 1;
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Called - size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" "<<geoStart<<","<<geoEnd ,10);
int maxGridSize = sizex; // get max size
if(sizey>maxGridSize) maxGridSize = sizey;
@@ -102,16 +119,19 @@ void initGridSizes(int &sizex, int &sizey, int &sizez,
for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }
// at least size 4 on coarsest level
int minSize = 2<<(maskBits+2);
int minSize = 4<<mMaxRefine; //(maskBits+2);
if(sizex<minSize) sizex = minSize;
if(sizey<minSize) sizey = minSize;
if(sizez<minSize) sizez = minSize;
sizeMask = ~sizeMask;
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask) ,10);
if(debugGridsizeInit) debMsgStd("initGridSizes",DM_MSG,"Size X:"<<sizex<<" Y:"<<sizey<<" Z:"<<sizez<<" m"<<convertFlags2String(sizeMask)<<", maskBits:"<<maskBits<<",minSize:"<<minSize ,10);
sizex &= sizeMask;
sizey &= sizeMask;
sizez &= sizeMask;
// debug
int nextinc = (~sizeMask)+1;
debMsgStd("initGridSizes",DM_MSG,"MPTeffMLtest, curr size:"<<PRINT_VEC(sizex,sizey,sizez)<<", next size:"<<PRINT_VEC(sizex+nextinc,sizey+nextinc,sizez+nextinc) ,10);
// force geom size to match rounded/modified grid sizes
geoEnd[0] = geoStart[0] + cellSize*(LbmFloat)sizex;
@@ -119,10 +139,11 @@ void initGridSizes(int &sizex, int &sizey, int &sizez,
geoEnd[2] = geoStart[2] + cellSize*(LbmFloat)sizez;
}
void calculateMemreqEstimate( int resx,int resy,int resz, int refine,
void calculateMemreqEstimate( int resx,int resy,int resz,
int refine, float farfield,
double *reqret, string *reqstr) {
// debug estimation?
const bool debugMemEst = false;
const bool debugMemEst = true;
// COMPRESSGRIDS define is not available here, make sure it matches
const bool useGridComp = true;
// make sure we can handle bid numbers here... all double
@@ -153,10 +174,17 @@ void calculateMemreqEstimate( int resx,int resy,int resz, int refine,
if(debugMemEst) debMsgStd("calculateMemreqEstimate",DM_MSG,"refine "<<i<<", mc:"<<memCnt, 10);
}
// isosurface memory
memCnt += (double)( (3*sizeof(int)+sizeof(float)) * ((resx+2)*(resy+2)*(resz+2)) );
// isosurface memory, use orig res values
if(farfield>0.) {
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:"<<memCnt, 10);
// cpdata init check missing...
double memd = memCnt;
char *sizeStr = "";
const double sfac = 1024.0;
@@ -187,7 +215,7 @@ void LbmSolverInterface::initDomainTrafo(float *mat) {
/*******************************************************************************/
/*! parse a boundary flag string */
CellFlagType LbmSolverInterface::readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed) {
string val = mpAttrs->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; i<mpGiObjects->size(); i++) {
if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 );
//debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"id:"<<mLbmInitId<<" obj:"<< (*mpGiObjects)[i]->getName() <<" gid:"<<(*mpGiObjects)[i]->getGeoInitId(), 9 );
}
debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Accurate geo init: "<<mAccurateGeoinit, 9)
debMsgStd("LbmSolverInterface::initGeoTree",DM_MSG,"Objects: "<<mpGiObjects->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<int> giObjFirstHistSide;
giObjFirstHistSide.resize( mpGiObjects->size() );
for(size_t i=0; i<mGiObjInside.size(); i++) {
@@ -335,7 +374,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) {
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);

61
intern/elbeem/intern/solver_interface.h
View File

@@ -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<ntlGeometryObject*> getDebugObjects() { vector<ntlGeometryObject*> 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 );

644
intern/elbeem/intern/solver_main.cpp
View File

@@ -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;j<mLevel[mMaxRefine].lSizey-0;j++) {
for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
//float val = 0.;
if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) {
//val = QCELL(mMaxRefine,i,j,k, workSet,dFfrac);
FORDF0{
const double df = (double)QCELL(mMaxRefine,i,j,k, workSet,l);
globdfavg[l] += df;
if(df>globdfmax[l]) globdfmax[l] = df;
//if(df<=0.01) { errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l); }
if(df<globdfmin[l]) globdfmin[l] = df;
//errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l<<" "<<df<<" "<<globdfmin[l]);
}
globdfcnt+=1.;
}
}
}
}
if(this->mStepCnt%10==0) {
FORDF0{ errMsg("GLOBDF","l="<<l<<" avg="<<(globdfavg[l]/globdfcnt)<<" max="<<globdfmax[l]<<" min="<<globdfmin[l]<<" "<<globdfcnt); }
}
}
}
// intlbm test */
stepMain();
}
void LbmFsgrSolver::stepMain()
{
this->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 : "<<this->mStepCnt, 10);
if(!this->mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime, 10); }
//debMsgDirect( "LbmFsgrSolver::step : "<<this->mStepCnt<<" ");
myTime_t timestart = getTime();
//debugOutNnl("LbmFsgrSolver::step : "<<mStepCnt, 10);
if(!mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime, 10); }
//debMsgDirect( "LbmFsgrSolver::step : "<<mStepCnt<<" ");
//myTime_t timestart = 0;
//if(mStartSymm) { checkSymmetry("step1"); } // DEBUG
@@ -93,15 +58,15 @@ void LbmFsgrSolver::stepMain()
// important - keep for tadap
LbmFloat lastMass = mCurrentMass;
mCurrentMass = this->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:"<<this->mStepCnt<<" s-1:"<<(this->mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
int dsbits = mStepCnt ^ (mStepCnt-1);
//errMsg("S"," step:"<<mStepCnt<<" s-1:"<<(mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
for(int lev=0; lev<=mMaxRefine; lev++) {
//if(! (this->mStepCnt&(1<<lev)) ) {
//if(! (mStepCnt&(1<<lev)) ) {
if( dsbits & (1<<(mMaxRefine-lev)) ) {
//errMsg("S"," l"<<lev);
@@ -124,50 +89,44 @@ void LbmFsgrSolver::stepMain()
}
// prepare next step
this->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:"<<this->mStepCnt<<" t:"<<mSimulationTime<<
//debMsgDirect(
" mlsups(curr:"<<this->mMLSUPS<<
" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
" totcls:"<<(this->mNumUsedCells)<< sepStr<<
" avgcls:"<< avgcls<< sepStr<<
" intd:"<<mNumInterdCells<< sepStr<<
" invif:"<<mNumInvIfCells<< sepStr<<
" invift:"<<mNumInvIfTotal<< sepStr<<
" fsgrcs:"<<mNumFsgrChanges<< sepStr<<
" filled:"<<this->mNumFilledCells<<", emptied:"<<this->mNumEmptiedCells<< sepStr<<
" mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
" probs:"<<mNumProblems<< sepStr<<
" simt:"<<mSimulationTime<< sepStr<<
" for '"<<this->mName<<"' " , 10);
if(!mSilent){
int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt);
debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime<<
" cur-mlsups:"<<mMLSUPS<< //" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<<
" totcls:"<<mNumUsedCells<< " avgcls:"<< avgcls<<
" intd:"<<mNumInterdCells<< " invif:"<<mNumInvIfCells<<
" invift:"<<mNumInvIfTotal<< " fsgrcs:"<<mNumFsgrChanges<<
" filled:"<<mNumFilledCells<<", emptied:"<<mNumEmptiedCells<<
" mMxv:"<<PRINT_VEC(mMxvx,mMxvy,mMxvz)<<", tscnts:"<<mTimeSwitchCounts<<
//" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
" probs:"<<mNumProblems<< " simt:"<<mSimulationTime<<
" took:"<< getTimeString(timeend-timestart)<<
" for '"<<mName<<"' " , 10);
} else { debMsgDirect("."); }
if(this->mStepCnt==1) {
mMinNoCells = mMaxNoCells = this->mNumUsedCells;
if(mStepCnt==1) {
mMinNoCells = mMaxNoCells = mNumUsedCells;
} else {
if(this->mNumUsedCells>mMaxNoCells) mMaxNoCells = this->mNumUsedCells;
if(this->mNumUsedCells<mMinNoCells) mMinNoCells = this->mNumUsedCells;
if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells;
if(mNumUsedCells<mMinNoCells) mMinNoCells = mNumUsedCells;
}
// mass scale test
@@ -191,7 +150,7 @@ void LbmFsgrSolver::stepMain()
errMsg("MDTDD","\n\n");
errMsg("MDTDD","FORCE RESCALE MASS! "
<<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
<<" step:"<<this->mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
<<" step:"<<mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
);
errMsg("MDTDD","\n\n");
@@ -233,8 +192,8 @@ void LbmFsgrSolver::stepMain()
}
#if LBM_INCLUDE_TESTSOLVERS==1
if((mUseTestdata)&&(this->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:"<<this->mStepCnt
<<": dccd="<< mCurrentMass
<<",d"<<mdelta
<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
<<" totst:"<<getTimeString(timeend-timestart), 3);
double effMLSUPS = ((double)mNumUsedCells / ((timeend2-timestart)/(double)1000.0) ) / (1000000.0);
if(mInitDone) {
if(effMLSUPS>10000){ effMLSUPS = -1; }
else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
}
debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<<glob_mpindex<<" step:"<<mStepCnt
<<" dccd="<< mCurrentMass
//<<",d"<<mdelta
//<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
//<<"/"<<mCurrentVolume<<"(fix="<<mFixMass<<",ini="<<mInitialMass<<"), "
<<" effMLSUPS=("<< effMLSUPS
<<",avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "
<<" took totst:"<< getTimeString(timeend2-timestart), 3);
// nicer output
debMsgDirect(std::endl);
//debMsgDirect(std::endl);
// */
messageOutputForce("");
//#endif // ELBEEM_PLUGIN!=1
}
@@ -291,15 +260,15 @@ void LbmFsgrSolver::fineAdvance()
// warning assume -Y gravity...
mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
if(mFVHeight<1.0) mFVHeight = 1.0;
this->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 "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
// update other set
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
@@ -308,7 +277,7 @@ void LbmFsgrSolver::fineAdvance()
// flag init... (work on current set, to simplify flag checks)
reinitFlags( mLevel[mMaxRefine].setCurr );
if(!this->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 "<<PRINT_IJK<<"="<<
RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
(*pFlagSrc)<<","<<(*pFlagDst) <<", pointers "<<
(int)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(int)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
(int)(pFlagSrc)<<","<<(int)(pFlagDst)<<" "
(long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
(long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
);
CAUSE_PANIC;
}
if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) ||
(&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
(int)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(int)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
(int)(ccel)<<","<<(int)(tcel)<<" "
(long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
(long)(ccel)<<","<<(long)(tcel)<<" "
);
CAUSE_PANIC;
}
@@ -466,7 +435,7 @@ LbmFsgrSolver::mainLoop(int lev)
// old INTCFCOARSETEST==1
if( (oldFlag & (CFGrFromCoarse)) ) {
if(( this->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( (i<cutMin)||(i>this->mSizex-cutMin)||
(j<cutMin)||(j>this->mSizey-cutMin)||
(k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
if( (i<cutMin)||(i>mSizex-cutMin)||
(j<cutMin)||(j>mSizey-cutMin)||
(k<cutMin)||(k>mSizez-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) && (lcsmqo<P_LCSMQO)) )// normal checks
||(k>this->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 "<<PRINT_IJK<<" zfac"<<zfac);
//if( (prob< (basethresh*rl)) && (lcsmqo>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:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
} // wind disturbance
// "wind" disturbance
// use realworld relative velocity here instead?
if( (doAdd &&
((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
||(k>mSizez-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 "<<PRINT_IJK<<" zfac"<<zfac);
} else {
// normal probability
//? LbmFloat fac = P_LCSMQO-lcsmqo;
//? jdf *= fac;
}
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:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
} // wind disturbance
} // mPartUsePhysModel
else {
// empirical model
//if((prob<basethresh) && (lcsmqo>0.0095)) { // add
if((prob<basethresh) && (lcsmqo>0.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 "<<PRINT_IJK<<"rwl:"<<rl<<" usqr:"<<usqr <<" qo:"<<lcsmqo<<" add="<<doAdd );
if( (doAdd) ) { // ADD DROP
LbmFloat len = norm(LbmVec(ux,uy,uz));
// WHMOD
//for(int s=0; s<10; s++) { // multiple parts
// dont try to subtract from empty cells
// ensure cell has enough mass for new drop
LbmFloat newPartsize = 1.0;
if(mPartUsePhysModel) {
// 1-10
newPartsize += 9.0* (rand()/(RAND_MAX+1.0));
} else {
// 1-5, overall size has to be less than
// .62 (ca. 0.5) to make drops significantly smaller
// than a full cell!
newPartsize += 4.0* (rand()/(RAND_MAX+1.0));
}
LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); //PARTMASS(mPartDropMassSub*newPartsize); // mass: 4/3 pi r^3 rho
while(dropmass>mass) {
newPartsize -= 0.2;
dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize);
}
if(newPartsize<=1.) doAdd=false;
if( (doAdd) ) { // init new particle
for(int s=0; s<1; s++) { // one part!
//LbmFloat prob = this->mPartGenProb * 0.02* (rand()/(RAND_MAX+1.0));
const LbmFloat posjitter = 1.0;
const LbmFloat posjitter = 0.05;
const LbmFloat posjitteroffs = posjitter*-0.5;
LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
const LbmFloat jitterstr = 0.1;
const LbmFloat jitterstr = 1.0;
const LbmFloat jitteroffs = jitterstr*-0.5;
LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
@@ -964,40 +949,46 @@ LbmFsgrSolver::mainLoop(int lev)
// average normal & velocity
// -> mostly along velocity dir, many into surface
LbmVec pv = (LbmVec(nx+jx,ny+jy,nz+jz)*0.75 + getNormalized(LbmVec(ux,uy,uz)) )*0.35;
normalize(pv);
// fluid velocity (not normalized!)
LbmVec flvelVel = LbmVec(ux,uy,uz);
LbmFloat flvelLen = norm(flvelVel);
// surface normal
LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]);
normalize(normVel);
LbmFloat normScale = (0.01+flvelLen);
// jitter vector, 0.2 * flvel
LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1;
// weighten velocities
const LbmFloat flvelWeight = 0.9;
LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel;
LbmFloat srci = i+0.5+jpx; // TEST? + (pv[0]*1.41);
LbmFloat srcj = j+0.5+jpy; // TEST? + (pv[1]*1.41);
LbmFloat srck = k+0.5+jpz; // TEST? + (pv[2]*1.41);
// offset towards surface (hide popping)
jpx += -normVel[0]*0.4;
jpy += -normVel[1]*0.4;
jpz += -normVel[2]*0.4;
LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz;
int type=0;
//if((s%3)!=2) {} else { type=PART_FLOAT; }
//type = PART_DROP;
type = PART_INTER;
// drop
/*srci += (pv[0]*1.41);
srcj += (pv[1]*1.41);
srck += (pv[2]*1.41);
if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction */
type = PART_DROP;
pv *= len;
LbmFloat size = 1.0+ 9.0* (rand()/(RAND_MAX+1.0));
mpParticles->addParticle(srci, srcj, srck); //i+0.5+jpx,j+0.5+jpy,k+0.5+jpz);
mpParticles->getLast()->setVel(pv[0],pv[1],pv[2]);
//? mpParticles->getLast()->advanceVel(); // advance a bit outwards
mpParticles->getLast()->setStatus(PART_IN);
mpParticles->getLast()->setType(type);
//if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
mpParticles->getLast()->setSize(size);
//errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
//mass -= size*0.001; // NTEST!
mass -= size*0.0020; // NTEST!
# if LBMDIM==2
mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
newpartVel[2]=0.; srck=0.5;
# endif // LBMDIM==2
//errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel() );
// 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 "<<PRINT_IJK<<" u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize );
//errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel );
FSGR_ADDPART(np);
} // new part
//errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel() );
} // multiple parts
} // doAdd
} // */
@@ -1047,7 +1038,7 @@ LbmFsgrSolver::mainLoop(int lev)
if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1; // NEWSURFT
filledp.x = i; filledp.y = j; filledp.z = k;
LIST_FULL(filledp);
//this->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 "<<PRINT_IJK<<" id"<<id);
GRID_REGION_START();
GRID_LOOP_START();
for(int l=0; l<dTotalNum; l++) { RAC(ccel,l) = 0.; }
*pFlagSrc =0;
*pFlagDst =0;
//errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
#if PARALLEL!=1
GRID_LOOPREG_END();
GRID_LOOPREG_END();
#else // PARALLEL==1
#include "paraloopend.h" // = GRID_LOOPREG_END();
#endif // PARALLEL==1
//GRID_REGION_END();
// TEST! */
/* dummy remove warnings */
calcCurrentMass = calcCurrentVolume = 0.;
calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
// change grid
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
mLevel[mMaxRefine].setCurr ^= 1;
/* dummy, remove warnings */
calcCurrentMass = calcCurrentVolume = 0.;
calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
// change grid
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
mLevel[mMaxRefine].setCurr ^= 1;
}
}
@@ -1176,14 +1158,9 @@ LbmFsgrSolver::standingFluidPreinit()
# endif // OPT3D==true
GRID_LOOP_START();
//FORDF0{ RAC(ccel,l) = 0.; }
//*pFlagSrc =0;
//*pFlagDst =0;
//errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
if( (currFlag & (CFEmpty|CFBnd)) ) continue;
//ccel = RACPNT(lev, i,j,k,workSet);
//tcel = RACPNT(lev, i,j,k,otherSet);
if( (currFlag & (CFInter)) ) {
// copy all values
@@ -1198,7 +1175,6 @@ LbmFsgrSolver::standingFluidPreinit()
nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
}
DEFAULT_STREAM;
//ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2];
DEFAULT_COLLIDEG(mLevel[lev].gravity);
}
for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
@@ -1207,8 +1183,6 @@ LbmFsgrSolver::standingFluidPreinit()
#else // PARALLEL==1
#include "paraloopend.h" // = GRID_LOOPREG_END();
#endif // PARALLEL==1
//GRID_REGION_END();
// TEST! */
/* dummy remove warnings */
calcCurrentMass = calcCurrentVolume = 0.;
@@ -1227,23 +1201,13 @@ LbmFsgrSolver::standingFluidPreinit()
LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
//return 0.0; // test
int level = mMaxRefine;
LbmFloat *ccel = RACPNT(level, i,j,k, workSet);
LbmFloat *ccel = RACPNT(level, i,j,k, workSet);
LbmFloat nx,ny,nz, nv1,nv2;
if(RFLAG_NB(level,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
if(RFLAG_NB(level,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
nx = 0.5* (nv2-nv1);
if(RFLAG_NB(level,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
if(RFLAG_NB(level,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
ny = 0.5* (nv2-nv1);
#if LBMDIM==3
if(RFLAG_NB(level,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
if(RFLAG_NB(level,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
nz = 0.5* (nv2-nv1);
#else //LBMDIM==3
nz = 0.0;
#endif //LBMDIM==3
LbmFloat scal = mDvecNrm[l][0]*nx + mDvecNrm[l][1]*ny + mDvecNrm[l][2]*nz;
// computenormal
CellFlagType *cflag = &RFLAG(level, i,j,k, workSet);
LbmFloat n[3];
computeFluidSurfaceNormal(ccel,cflag, n);
LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2];
LbmFloat ret = 1.0;
// forward direction, add mass (for filling cells):
@@ -1259,6 +1223,115 @@ LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workS
return ret;
}
// warning - normal compuations are without
// boundary checks &
// normalization
void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) {
const int level = mMaxRefine;
LbmFloat nx,ny,nz, nv1,nv2;
const CellFlagType flagE = *(cflagpnt+1);
const CellFlagType flagW = *(cflagpnt-1);
if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); }
else if(flagE &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); }
else if(flagW &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
nx = 0.5* (nv2-nv1);
const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); }
else if(flagN &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); }
else if(flagS &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
ny = 0.5* (nv2-nv1);
#if LBMDIM==3
const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); }
else if(flagT &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); }
else if(flagB &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
nz = 0.5* (nv2-nv1);
#else //LBMDIM==3
nz = 0.0;
#endif //LBMDIM==3
// return vals
snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
LbmFloat nx=0.,ny=0.,nz=0.;
ccel = NULL; cflagpnt=NULL; // remove warning
snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
const int level = mMaxRefine;
LbmFloat nx,ny,nz, nv1,nv2;
ccel = NULL; // remove warning
const CellFlagType flagE = *(cflagpnt+1);
const CellFlagType flagW = *(cflagpnt-1);
if(flagE &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagW &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
nx = 0.5* (nv2-nv1);
const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
if(flagN &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagS &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
ny = 0.5* (nv2-nv1);
#if LBMDIM==3
const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
if(flagT &(CFBnd)){ nv1 = 1.; }
else nv1 = 0.0;
if(flagB &(CFBnd)){ nv2 = 1.; }
else nv2 = 0.0;
nz = 0.5* (nv2-nv1);
#else //LBMDIM==3
nz = 0.0;
#endif //LBMDIM==3
// return vals
snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) {
bool nonorm = false;
LbmFloat nx=0.,ny=0.,nz=0.;
if(i<=0) { nx = 1.; nonorm = true; }
if(i>=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<<" "<<this->mFixMass);
//errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass);
}
weightIndex++;
@@ -1531,9 +1604,9 @@ void LbmFsgrSolver::reinitFlags( int workSet ) {
massChange = 0.0;
} else {
// Problem! no interface neighbors
this->mFixMass += massChange;
mFixMass += massChange;
//TTT mNumProblems++;
//errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<this->mFixMass);
//errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass);
}
weightIndex++;
@@ -1582,7 +1655,7 @@ void LbmFsgrSolver::reinitFlags( int workSet ) {
continue;
} // */
QCELL(workLev,i,j,k, workSet, dMass) += (this->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:"<<this->mFixMass<<" nif:"<<mListNewInter.size() );
this->mFixMass = 0.0;
//errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() );
mFixMass = 0.0;
}
// empty lists for next step

13
intern/elbeem/intern/solver_relax.h
View File

@@ -3,7 +3,7 @@
* El'Beem - the visual lattice boltzmann freesurface simulator
* 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
*
* Combined 2D/3D Lattice Boltzmann relaxation macros
*
@@ -378,6 +378,15 @@
} \
} else { \
m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
if(RFLAG(lev, i,j,k, mLevel[lev].setCurr)&CFFluid) { \
if(!(nbf&(CFFluid|CFInter)) ) { \
int ni=i+this->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"<<l<<" value:"<<m[l]<<" at "<<PRINT_IJK<<" from "<< \
PRINT_VEC(ni,nj,nk) <<" l"<<(l)<< \
" nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther) ); \
\
\
} } \
STREAMCHECK(4, i+this->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
};

774
intern/elbeem/intern/solver_util.cpp
View File

File diff suppressed because it is too large Load Diff

154
intern/elbeem/intern/utilities.cpp
View File

@@ -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 <png.h>
#endif
#endif // NOPNG
#include <zlib.h>
// 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;j<h;j++) {
for(int i=0;i<h;i++) {
// remove alpha values
gzwrite(gzf,&rowsp[j][i*4],3);
}
}
gzclose( gzf );
errMsg("writePng/ppm","Write_png/ppm: wrote to "<<filennew.str()<<".");
return 0;
fail:
errMsg("writePng/ppm","Write_png/ppm: could not write to "<<filennew.str()<<" !");
return -1;
}
#endif // NOPNG
@@ -243,7 +290,21 @@ static string col_purple ( "\033[0;35m");
static string col_bright_purple ( "\033[1;35m");
static string col_neutral ( "\033[0m");
static string col_std = col_bright_gray;
std::ostringstream globOutstr;
bool globOutstrForce=false;
#define DM_NONE 100
void messageOutputForce(string from) {
bool org = globOutstrForce;
globOutstrForce = true;
messageOutputFunc(from, DM_NONE, "\n", 0);
globOutstrForce = org;
}
void messageOutputFunc(string from, int id, string msg, myTime_t interval) {
// fast skip
if((id!=DM_FATAL)&&(gDebugLevel<=0)) return;
if(interval>0) {
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("<<gElbeemState<<"):" << col_red;
break;
case DM_NONE:
// only internal debugging msgs
break;
default:
// this shouldnt happen...
sout << col_red << " --- messageOutputFunc error: invalid id ("<<id<<") --- aborting... \n\n" << col_std;
@@ -303,19 +369,61 @@ void messageOutputFunc(string from, int id, string msg, myTime_t interval) {
sout << "\n"; // add newline for output
}
#ifdef WIN32
// debug level is >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 '"<<strEnvName<<"'='"<<getenv(strEnvName)<<"', debugLevel set to: "<<gDebugLevel<<"\n", 1);
}
if(getenv(strEnvName2)) {
gDebugLevel = atoi(getenv(strEnvName2));
if(gDebugLevel>0) debMsgStd("performElbeemSimulation",DM_NOTIFY,"Using envvar '"<<strEnvName2<<"'='"<<getenv(strEnvName2)<<"', debugLevel set to: "<<gDebugLevel<<"\n", 1);
}
if(gDebugLevel< 0) gDebugLevel = 0;
if(gDebugLevel>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

6
intern/elbeem/intern/utilities.h
View File

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

4
source/blender/blenkernel/intern/DerivedMesh.c
View File

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

14
source/blender/blenkernel/intern/effect.c
View File

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

3
source/blender/makesdna/DNA_object_fluidsim.h
View File

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

82
source/blender/src/buttons_object.c
View File

@@ -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)<MAX_EFFECT) &&
(!paf)) {
// create new entry
copy_act_effect(ob); DAG_scene_sort(G.scene); }
paf = give_parteff(ob);
paf->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 {

296
source/blender/src/fluidsim.c
View File

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