- bugfixes
authorNils Thuerey <nils@thuerey.de>
Sun, 5 Nov 2006 16:30:29 +0000 (16:30 +0000)
committerNils Thuerey <nils@thuerey.de>
Sun, 5 Nov 2006 16:30:29 +0000 (16:30 +0000)
  #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

48 files changed:
intern/elbeem/COPYING
intern/elbeem/extern/elbeem.h
intern/elbeem/intern/attributes.cpp
intern/elbeem/intern/attributes.h
intern/elbeem/intern/elbeem.cpp
intern/elbeem/intern/elbeem.h
intern/elbeem/intern/isosurface.cpp
intern/elbeem/intern/isosurface.h
intern/elbeem/intern/loop_tools.h
intern/elbeem/intern/ntl_blenderdumper.cpp
intern/elbeem/intern/ntl_blenderdumper.h
intern/elbeem/intern/ntl_bsptree.cpp
intern/elbeem/intern/ntl_bsptree.h
intern/elbeem/intern/ntl_geometryclass.h
intern/elbeem/intern/ntl_geometrymodel.cpp
intern/elbeem/intern/ntl_geometrymodel.h
intern/elbeem/intern/ntl_geometryobject.cpp
intern/elbeem/intern/ntl_geometryobject.h
intern/elbeem/intern/ntl_geometryshader.h
intern/elbeem/intern/ntl_lighting.cpp
intern/elbeem/intern/ntl_lighting.h
intern/elbeem/intern/ntl_matrices.h
intern/elbeem/intern/ntl_ray.cpp
intern/elbeem/intern/ntl_ray.h
intern/elbeem/intern/ntl_vector3dim.h
intern/elbeem/intern/ntl_world.cpp
intern/elbeem/intern/ntl_world.h
intern/elbeem/intern/parametrizer.cpp
intern/elbeem/intern/parametrizer.h
intern/elbeem/intern/particletracer.cpp
intern/elbeem/intern/particletracer.h
intern/elbeem/intern/simulation_object.cpp
intern/elbeem/intern/simulation_object.h
intern/elbeem/intern/solver_adap.cpp
intern/elbeem/intern/solver_class.h
intern/elbeem/intern/solver_init.cpp
intern/elbeem/intern/solver_interface.cpp
intern/elbeem/intern/solver_interface.h
intern/elbeem/intern/solver_main.cpp
intern/elbeem/intern/solver_relax.h
intern/elbeem/intern/solver_util.cpp
intern/elbeem/intern/utilities.cpp
intern/elbeem/intern/utilities.h
source/blender/blenkernel/intern/DerivedMesh.c
source/blender/blenkernel/intern/effect.c
source/blender/makesdna/DNA_object_fluidsim.h
source/blender/src/buttons_object.c
source/blender/src/fluidsim.c

index d876362..2600c73 100644 (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
index ea21e80..b3feda8 100644 (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
index b7f1945..0cf5315 100644 (file)
 /******************************************************************************
  *
  * 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);
-}
-
-//! get channel as integer value
+       return AnimChannel<double>();
+}
 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);
-}
-
-//! get channel as integer value
+       return AnimChannel<int>();
+}
 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);
-}
-
-//! get channel as float vector set
+       return AnimChannel<ntlVec3d>();
+}
 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;
index 20e5341..1935509 100644 (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,178 +110,84 @@ class AnimChannel
 };
 
 
-//! A single attribute
+// helper class (not templated) for animated meshes
+class ntlSetVec3f {
+       public:
+               ntlSetVec3f(): mVerts() {};
+               ntlSetVec3f(double v);
+               ntlSetVec3f(vector<ntlVec3f> &v) { mVerts = v; };
+
+               const ntlSetVec3f& operator=(double v );
+               ntlSetVec3f& operator+=( double v );
+               ntlSetVec3f& operator+=( const ntlSetVec3f &v );
+               ntlSetVec3f& operator*=( double v );
+               ntlSetVec3f& operator*=( const ntlSetVec3f &v );
+
+               vector<ntlVec3f> mVerts;
+};
+
+
+// warning: DEPRECATED - replaced by elbeem API
 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 
+       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);
-               //! 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:
-               ntlSetVec3f(): mVerts() {};
-               ntlSetVec3f(double v);
-               ntlSetVec3f(vector<ntlVec3f> &v) { mVerts = v; };
-
-               const ntlSetVec3f& operator=(double v );
-               ntlSetVec3f& operator+=( double v );
-               ntlSetVec3f& operator+=( const ntlSetVec3f &v );
-               ntlSetVec3f& operator*=( double v );
-               ntlSetVec3f& operator*=( const ntlSetVec3f &v );
-
-               vector<ntlVec3f> mVerts;
-};
-
+// 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
 
index d4242da..b2779f5 100644 (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]) );
index ea21e80..b3feda8 100644 (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
index 283fd17..fdef04d 100644 (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);
-
-                               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);
-
-                               /*int bndskip = 0; // BNDOFFT
-                               for(int s=0; s<8; s++) if(value[s]==-76.) bndskip++;
-                               if(bndskip>0) continue; // */
-
-                               // 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;
-                               }
+       if(mSubdivs<=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;
+
+                                       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;
+
+                                       // 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);
+                                       // 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 ];
+
+                                       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 
+                                       }
 
-                                                       //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 { 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!");
+               }
 
-                               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;
+               // 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() );
+                                                       }
 
-                               // 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] ] );
+                                                       } // 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;
                                }
-                               
-      }//i
-    }// j
-  } // k
-
-  // 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 );
+                       // */
+
+               } // 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?
-
-               //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];
-               // 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];
+               if(d2 >= 9.0f) continue;
+
+               // aggressive smoothing factor
+               float smstr = 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
 }
 
 
index 5b51981..9902b40 100644 (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:
 
index 927263b..3c15a60 100644 (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 */ \
index 1bdea59..b1fece2 100644 (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 {
index 57f155f..df66ad7 100644 (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
  *
index afa1fa5..743ea79 100644 (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
  *
index face7c3..35bc7c6 100644 (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
  *
index 0b7a4e9..9282371 100644 (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 */
index ba6e638..0d769cf 100644 (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
  *
index 95cadc4..572440f 100644 (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
  *
index d815fb5..bc004b6 100644 (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;
 }
index a9a8109..a5afd6b 100644 (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:
 
index aba6b59..3ecb82e 100644 (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
  *
index d4eae70..cc0e6db 100644 (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
  *
index af53a13..772f01e 100644 (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
index 7d27a6e..1dd6159 100644 (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
index c838e6f..a25f420 100644 (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
  *
index 8c0188a..226c5f5 100644 (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
  *
index ba3ea68..d6a7557 100644 (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
  *
index befcbed..835ca7d 100644 (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(),
index 18f388a..9c5324c 100644 (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
  *
index b9dc1c3..ae2bc7f 100644 (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);
index c0df399..e05db12 100644 (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;
index 1d15cec..c537a89 100644 (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
index f08bd58..aae92aa 100644 (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();
index 3f954f1..fcec398 100644 (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++;
                }
+               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++) {
-                       if( flag&(1<<j) ) flagCount[j]++;
+                       std::ostringstream out;
+                       if(flagCount[j]>0) {
+                               out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1<<j) ) <<", " << charNl;
+                               debugOutNnl(out.str(), 5);
+                       }
                }
-               totalCells++;
-       }
-       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;
+               // 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...
index 200cad8..56d7f20 100644 (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
  *
index 9989208..ef516a5 100644 (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);
+
+       // 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
 
-       //newdt = this->mpParam->getTimestep() * (allowMax/nextmax);
-       LbmFloat newdt = this->mpParam->getTimestep(); // newtr
+       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");
                        }
index 578e107..e6b1ad4 100644 (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
index 1a19e1d..6f89131 100644 (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
  *
@@ -13,9 +13,6 @@
 // for geo init FGI_ defines
 #include "elbeem.h"
 
-// all boundary types at once
-#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
-
 // helper for 2d init
 #define SWAPYZ(vec) { \
                const LbmFloat tmp = (vec)[2]; \
        };
 
        /* precalculated equilibrium dfs, inited in lbmsolver constructor */
-       LbmFloat LbmFsgrSolver::dfEquil[ cDfNum ];
+       LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
 
 #else // end LBMDIM==3 , LBMDIM==2
 
        };
 
        /* precalculated equilibrium dfs, inited in lbmsolver constructor */
-       LbmFloat LbmFsgrSolver::dfEquil[ cDfNum ];
+       LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
 
 // D2Q9 end
 #endif  // LBMDIM==2
 
 
+// required globals
+extern bool glob_mpactive;
+extern int glob_mpnum, glob_mpindex;
+
+
 /******************************************************************************
  * Lbm Constructor
  *****************************************************************************/
@@ -318,7 +320,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mIsoWeightMethod(1),
        mMaxRefine(1), 
        mDfScaleUp(-1.0), mDfScaleDown(-1.0),
-       mInitialCsmago(0.03), // set to 0.02 for mMaxRefine==0 below
+       mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03
        mDebugOmegaRet(0.0),
        mLastOmega(1e10), mLastGravity(1e10),
        mNumInvIfTotal(0), mNumFsgrChanges(0),
@@ -329,38 +331,45 @@ LbmFsgrSolver::LbmFsgrSolver() :
   // not much to do here... 
 #if LBM_INCLUDE_TESTSOLVERS==1
        mpTest = new LbmTestdata();
-#endif // ELBEEM_PLUGIN!=1
-       this->mpIso = new IsoSurface( this->mIsoValue );
+       mMpNum = mMpIndex = 0;
+       mOrgSizeX = 0;
+       mOrgStartX = 0.;
+       mOrgEndX = 0.;
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
+       mpIso = new IsoSurface( mIsoValue );
 
   // init equilibrium dist. func 
   LbmFloat rho=1.0;
   FORDF0 {
-               this->dfEquil[l] = this->getCollideEq( l,rho,  0.0, 0.0, 0.0);
+               dfEquil[l] = this->getCollideEq( l,rho,  0.0, 0.0, 0.0);
   }
+       dfEquil[dMass] = 1.;
+       dfEquil[dFfrac] = 1.;
+       dfEquil[dFlux] = FLUX_INIT;
 
        // init LES
        int odm = 0;
        for(int m=0; m<LBMDIM; m++) { 
-               for(int l=0; l<this->cDfNum; l++) { 
+               for(int l=0; l<cDfNum; l++) { 
                        this->lesCoeffDiag[m][l] = 
                        this->lesCoeffOffdiag[m][l] = 0.0;
                }
        }
        for(int m=0; m<LBMDIM; m++) { 
                for(int n=0; n<LBMDIM; n++) { 
-                       for(int l=1; l<this->cDfNum; l++) { 
+                       for(int l=1; l<cDfNum; l++) { 
                                LbmFloat em;
                                switch(m) {
-                                       case 0: em = this->dfDvecX[l]; break;
-                                       case 1: em = this->dfDvecY[l]; break;
-                                       case 2: em = this->dfDvecZ[l]; break;
+                                       case 0: em = dfDvecX[l]; break;
+                                       case 1: em = dfDvecY[l]; break;
+                                       case 2: em = dfDvecZ[l]; break;
                                        default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
                                }
                                LbmFloat en;
                                switch(n) {
-                                       case 0: en = this->dfDvecX[l]; break;
-                                       case 1: en = this->dfDvecY[l]; break;
-                                       case 2: en = this->dfDvecZ[l]; break;
+                                       case 0: en = dfDvecX[l]; break;
+                                       case 1: en = dfDvecY[l]; break;
+                                       case 2: en = dfDvecZ[l]; break;
                                        default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
                                }
                                const LbmFloat coeff = em*en;
@@ -383,7 +392,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mDvecNrm[0] = LbmVec(0.0);
   FORDF1 {
                mDvecNrm[l] = getNormalized( 
-                       LbmVec(this->dfDvecX[this->dfInv[l]], this->dfDvecY[this->dfInv[l]], this->dfDvecZ[this->dfInv[l]] ) 
+                       LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] ) 
                        ) * -1.0; 
        }
 
@@ -395,15 +404,15 @@ LbmFsgrSolver::LbmFsgrSolver() :
 #if ELBEEM_PLUGIN!=1
        errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
 #endif
-       for(int n=0;(n<this->cDirNum); n++) { mGaussw[n] = 0.0; }
-       //for(int n=0;(n<this->cDirNum); n++) { 
-       for(int n=0;(n<this->cDfNum); n++) { 
-               const LbmFloat d = norm(LbmVec(this->dfVecX[n], this->dfVecY[n], this->dfVecZ[n]));
+       for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
+       //for(int n=0;(n<cDirNum); n++) { 
+       for(int n=0;(n<cDfNum); n++) { 
+               const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
                LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
                mGaussw[n] = w;
                totGaussw += w;
        }
-       for(int n=0;(n<this->cDirNum); n++) { 
+       for(int n=0;(n<cDirNum); n++) { 
                mGaussw[n] = mGaussw[n]/totGaussw;
        }
 
@@ -414,7 +423,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
 /*****************************************************************************/
 LbmFsgrSolver::~LbmFsgrSolver()
 {
-  if(!this->mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; }
+  if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; }
 #if COMPRESSGRIDS==1
        delete [] mLevel[mMaxRefine].mprsCells[1];
        mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
@@ -426,16 +435,13 @@ LbmFsgrSolver::~LbmFsgrSolver()
                        if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
                }
        }
-       delete this->mpIso;
+       delete mpIso;
        if(mpPreviewSurface) delete mpPreviewSurface;
-
-#if LBM_INCLUDE_TESTSOLVERS==1
        // cleanup done during scene deletion...
-#endif // ELBEEM_PLUGIN!=1
 
        // always output performance estimate
        debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
-  if(!this->mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
+  if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
 }
 
 
@@ -449,19 +455,19 @@ void LbmFsgrSolver::parseAttrList()
        LbmSolverInterface::parseStdAttrList();
 
        string matIso("default");
-       matIso = this->mpAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
-       this->mpIso->setMaterialName( matIso );
-       this->mOutputSurfacePreview = this->mpAttrs->readInt("surfacepreview", this->mOutputSurfacePreview, "SimulationLbm","this->mOutputSurfacePreview", false );
-       mTimeAdap = this->mpAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
-       this->mDomainBound = this->mpAttrs->readString("domainbound", this->mDomainBound, "SimulationLbm","mDomainBound", false );
-       this->mDomainPartSlipValue = this->mpAttrs->readFloat("domainpartslip", this->mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );
-
-       mIsoWeightMethod= this->mpAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
-       mInitSurfaceSmoothing = this->mpAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
-       this->mSmoothSurface = this->mpAttrs->readFloat("smoothsurface", this->mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
-       this->mSmoothNormals = this->mpAttrs->readFloat("smoothnormals", this->mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
-
-       mFsSurfGenSetting = this->mpAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
+       matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
+       mpIso->setMaterialName( matIso );
+       mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
+       mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
+       mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false );
+       mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );
+
+       mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
+       mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
+       mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
+       mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
+
+       mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
        if(mFsSurfGenSetting==-1) {
                // all on
                mFsSurfGenSetting = 
@@ -470,46 +476,61 @@ void LbmFsgrSolver::parseAttrList()
        }
 
        // refinement
-       mMaxRefine = this->mRefinementDesired;
-       mMaxRefine  = this->mpAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
+       mMaxRefine = mRefinementDesired;
+       mMaxRefine  = mpSifAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
        if(mMaxRefine<0) mMaxRefine=0;
        if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
-       mDisableStandingFluidInit = this->mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
-       mInit2dYZ = this->mpAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
-       mForceTadapRefine = this->mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
+       mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
+       mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
+       mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
 
        // demo mode settings
-       mFVHeight = this->mpAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
+       mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
        // FIXME check needed?
-       mFVArea   = this->mpAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
+       mFVArea   = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
 
        // debugging - skip some time...
        double starttimeskip = 0.;
-       starttimeskip = this->mpAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
+       starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
        mSimulationTime += starttimeskip;
        if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
 
 #if LBM_INCLUDE_TESTSOLVERS==1
        mUseTestdata = 0;
-       mUseTestdata = this->mpAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
-       mpTest->parseTestdataAttrList(this->mpAttrs);
+       mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
+       mpTest->parseTestdataAttrList(mpSifAttrs);
 #ifdef ELBEEM_PLUGIN
        mUseTestdata=1; // DEBUG
 #endif // ELBEEM_PLUGIN
-       errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
 
+       mMpNum = mpSifAttrs->readInt("mpnum",  mMpNum ,"LbmFsgrSolver", "mMpNum", false);
+       mMpIndex = mpSifAttrs->readInt("mpindex",  mMpIndex ,"LbmFsgrSolver", "mMpIndex", false);
+       if(glob_mpactive) {
+               // used instead...
+               mMpNum = glob_mpnum;
+               mMpIndex = glob_mpindex;
+       } else {
+               glob_mpnum = mMpNum;
+               glob_mpindex = 0;
+       }
+       errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<<glob_mpactive<<", "<<glob_mpindex<<"/"<<glob_mpnum);
+       if(mMpNum>0) {
+               mUseTestdata=1; // needed in this case...
+       }
+
+       errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
 #else // LBM_INCLUDE_TESTSOLVERS!=1
-       // off by default
+       // not testsolvers, off by default
        mUseTestdata = 0;
        if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
 #endif // LBM_INCLUDE_TESTSOLVERS!=1
   //if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
 
        if(mMaxRefine==0) mInitialCsmago=0.02;
-       mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
+       mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
        // deprecated!
        float mInitialCsmagoCoarse = 0.0;
-       mInitialCsmagoCoarse = this->mpAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
+       mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
 #if USE_LES==1
 #else // USE_LES==1
        debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
@@ -524,15 +545,15 @@ void LbmFsgrSolver::parseAttrList()
 void LbmFsgrSolver::initLevelOmegas()
 {
        // no explicit settings
-       this->mOmega = this->mpParam->calculateOmega(mSimulationTime);
-       this->mGravity = vec2L( this->mpParam->calculateGravity(mSimulationTime) );
-       this->mSurfaceTension = 0.; //this->mpParam->calculateSurfaceTension(); // unused
+       mOmega = mpParam->calculateOmega(mSimulationTime);
+       mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) );
+       mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused
        if(mInit2dYZ) { SWAPYZ(mGravity); }
 
        // check if last init was ok
-       LbmFloat gravDelta = norm(this->mGravity-mLastGravity);
-       //errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<this->mOmega<<" - lom:"<<mLastOmega<<" gv:"<<this->mGravity<<" - "<<mLastGravity<<" , "<<gravDelta  );
-       if((this->mOmega == mLastOmega) && (gravDelta<=0.0)) return;
+       LbmFloat gravDelta = norm(mGravity-mLastGravity);
+       //errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<mOmega<<" - lom:"<<mLastOmega<<" gv:"<<mGravity<<" - "<<mLastGravity<<" , "<<gravDelta  );
+       if((mOmega == mLastOmega) && (gravDelta<=0.0)) return;
 
        if(mInitialCsmago<=0.0) {
                if(OPT3D==1) {
@@ -544,47 +565,49 @@ void LbmFsgrSolver::initLevelOmegas()
        // use Tau instead of Omega for calculations
        { // init base level
                int i = mMaxRefine;
-               mLevel[i].omega    = this->mOmega;
-               mLevel[i].timestep = this->mpParam->getTimestep();
+               mLevel[i].omega    = mOmega;
+               mLevel[i].timestep = mpParam->getTimestep();
                mLevel[i].lcsmago = mInitialCsmago; //CSMAGO_INITIAL;
                mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
                mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
        }
 
        // init all sub levels
+       LbmFloat coarseCsmago = mInitialCsmago;
+       coarseCsmago = 0.04; // try stabilizing
        for(int i=mMaxRefine-1; i>=0; i--) {
                //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
                double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
                nomega                = 1.0/nomega;
                mLevel[i].omega       = (LbmFloat)nomega;
                mLevel[i].timestep    = 2.0 * mLevel[i+1].timestep;
-               mLevel[i].lcsmago = mInitialCsmago;
+               mLevel[i].lcsmago = coarseCsmago;
                mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
                mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
        }
        
        // for lbgk
-       mLevel[ mMaxRefine ].gravity = this->mGravity / mLevel[ mMaxRefine ].omega;
+       mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
        for(int i=mMaxRefine-1; i>=0; i--) {
                // should be the same on all levels...
                // for lbgk
                mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
        }
 
-       mLastOmega = this->mOmega;
-       mLastGravity = this->mGravity;
+       mLastOmega = mOmega;
+       mLastGravity = mGravity;
        // debug? invalidate old values...
-       this->mGravity = -100.0;
-       this->mOmega = -100.0;
+       mGravity = -100.0;
+       mOmega = -100.0;
 
        for(int i=0; i<=mMaxRefine; i++) {
-               if(!this->mSilent) {
+               if(!mSilent) {
                        errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz 
                                        <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
                                        <<" cmsagp:"<<mLevel[i].lcsmago<<", "
                                        << " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
                } else {
-                       if(!this->mInitDone) {
+                       if(!mInitDone) {
                                debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
                                                <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
                        }
@@ -604,55 +627,75 @@ void LbmFsgrSolver::initLevelOmegas()
 /*! finish the init with config file values (allocate arrays...) */
 bool LbmFsgrSolver::initializeSolverMemory()
 {
-  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<this->mInitDone<<" "<<(void*)this,1);
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<mInitDone<<" "<<(void*)this,1);
 
        // init cppf stage
        if(mCppfStage>0) {
-               this->mSizex *= mCppfStage;
-               this->mSizey *= mCppfStage;
-               this->mSizez *= mCppfStage;
+               mSizex *= mCppfStage;
+               mSizey *= mCppfStage;
+               mSizez *= mCppfStage;
        }
 
        // size inits to force cubic cells and mult4 level dimensions
        // and make sure we dont allocate too much...
        bool memOk = false;
-       int orgSx = this->mSizex;
-       int orgSy = this->mSizey;
-       int orgSz = this->mSizez;
+       int orgSx = mSizex;
+       int orgSy = mSizey;
+       int orgSz = mSizez;
        double sizeReduction = 1.0;
        double memEstFromFunc = -1.0;
        string memreqStr("");   
+       bool firstMInit = true;
+       int minitTries=0;
        while(!memOk) {
-               initGridSizes( this->mSizex, this->mSizey, this->mSizez,
-                               this->mvGeoStart, this->mvGeoEnd, mMaxRefine, PARALLEL);
-               calculateMemreqEstimate( this->mSizex, this->mSizey, this->mSizez, mMaxRefine, &memEstFromFunc, &memreqStr );
+               minitTries++;
+               initGridSizes( mSizex, mSizey, mSizez,
+                               mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL);
+
+               // MPT
+#if LBM_INCLUDE_TESTSOLVERS==1
+               if(firstMInit) {
+                       mrSetup();
+               }
+#endif // LBM_INCLUDE_TESTSOLVERS==1
+               firstMInit=false;
+
+               calculateMemreqEstimate( mSizex, mSizey, mSizez, 
+                               mMaxRefine, mFarFieldSize, &memEstFromFunc, &memreqStr );
                
                double memLimit;
+               string memLimStr("-");
                if(sizeof(void*)==4) {
                        // 32bit system, limit to 2GB
                        memLimit = 2.0* 1024.0*1024.0*1024.0;
+                       memLimStr = string("2GB");
                } else {
                        // 64bit, just take 16GB as limit for now...
                        memLimit = 16.0* 1024.0*1024.0*1024.0;
+                       memLimStr = string("16GB");
                }
                if(memEstFromFunc>memLimit) {
                        sizeReduction *= 0.9;
-                       this->mSizex = (int)(orgSx * sizeReduction);
-                       this->mSizey = (int)(orgSy * sizeReduction);
-                       this->mSizez = (int)(orgSz * sizeReduction);
-                       debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<memEstFromFunc<<"/"<<memLimit<<", retrying: "
-                                       <<this->mSizex<<" Y:"<<this->mSizey<<" Z:"<<this->mSizez, 3 );
+                       mSizex = (int)(orgSx * sizeReduction);
+                       mSizey = (int)(orgSy * sizeReduction);
+                       mSizez = (int)(orgSz * sizeReduction);
+                       debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
+                                       //memEstFromFunc<<"/"<<memLimit<<", "<<
+                                       memreqStr<<"/"<<memLimStr<<", "<<
+                                       "retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
+                                       , 3 );
                } else {
                        memOk = true;
                } 
        }
        
-       this->mPreviewFactor = (LbmFloat)this->mOutputSurfacePreview / (LbmFloat)this->mSizex;
-  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<this->mSizex<<" Y:"<<this->mSizey<<" Z:"<<this->mSizez<<
-         ", Domain: "<<this->mvGeoStart<<":"<<this->mvGeoEnd<<", "<<(this->mvGeoEnd-this->mvGeoStart)<<
-         ", Intbytes: "<< sizeof(void*) <<
+       mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
+         ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
+         ", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
          ", est. Mem.Req.: "<<memreqStr        ,2);
-       this->mpParam->setSize(this->mSizex, this->mSizey, this->mSizez);
+       mpParam->setSize(mSizex, mSizey, mSizez);
+       if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }
 
   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
                <<"LBM_EPSILON="<<LBM_EPSILON  <<" "
@@ -666,22 +709,20 @@ bool LbmFsgrSolver::initializeSolverMemory()
                <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" " 
                <<"USE_LES="<<USE_LES <<" " 
                ,10);
-#if ELBEEM_PLUGIN!=1
-#endif // ELBEEM_PLUGIN!=1
 
        // perform 2D corrections...
-       if(LBMDIM == 2) this->mSizez = 1;
+       if(LBMDIM == 2) mSizez = 1;
 
-       this->mpParam->setSimulationMaxSpeed(0.0);
-       if(mFVHeight>0.0) this->mpParam->setFluidVolumeHeight(mFVHeight);
-       this->mpParam->setTadapLevels( mMaxRefine+1 );
+       mpParam->setSimulationMaxSpeed(0.0);
+       if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
+       mpParam->setTadapLevels( mMaxRefine+1 );
 
        if(mForceTadapRefine>mMaxRefine) {
-               this->mpParam->setTadapLevels( mForceTadapRefine+1 );
+               mpParam->setTadapLevels( mForceTadapRefine+1 );
                debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
        }
 
-       if(!this->mpParam->calculateAllMissingValues(mSimulationTime, false)) {
+       if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
                errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
                return false;
        }
@@ -706,9 +747,9 @@ bool LbmFsgrSolver::initializeSolverMemory()
        }
 
        // init sizes
-       mLevel[mMaxRefine].lSizex = this->mSizex;
-       mLevel[mMaxRefine].lSizey = this->mSizey;
-       mLevel[mMaxRefine].lSizez = this->mSizez;
+       mLevel[mMaxRefine].lSizex = mSizex;
+       mLevel[mMaxRefine].lSizey = mSizey;
+       mLevel[mMaxRefine].lSizez = mSizez;
        for(int i=mMaxRefine-1; i>=0; i--) {
                mLevel[i].lSizex = mLevel[i+1].lSizex/2;
                mLevel[i].lSizey = mLevel[i+1].lSizey/2;
@@ -722,8 +763,8 @@ bool LbmFsgrSolver::initializeSolverMemory()
        }
 
        double ownMemCheck = 0.0;
-       mLevel[ mMaxRefine ].nodeSize = ((this->mvGeoEnd[0]-this->mvGeoStart[0]) / (LbmFloat)(this->mSizex));
-       mLevel[ mMaxRefine ].simCellSize = this->mpParam->getCellSize();
+       mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
+       mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
        mLevel[ mMaxRefine ].lcellfactor = 1.0;
        LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
        // +4 for safety ?
@@ -759,8 +800,14 @@ bool LbmFsgrSolver::initializeSolverMemory()
                ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
        }
 
-       // isosurface memory
-       ownMemCheck += (3*sizeof(int)+sizeof(float)) * ((this->mSizex+2)*(this->mSizey+2)*(this->mSizez+2));
+       // isosurface memory, use orig res values
+       if(mFarFieldSize>0.) {
+               ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
+       } else {
+               // ignore 3 int slices...
+               ownMemCheck += (double)( (              sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
+       }
+
        // sanity check
 #if ELBEEM_PLUGIN!=1
        if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
@@ -782,31 +829,31 @@ bool LbmFsgrSolver::initializeSolverMemory()
 
        // calc omega, force for all levels
        initLevelOmegas();
-       mMinTimestep = this->mpParam->getTimestep();
-       mMaxTimestep = this->mpParam->getTimestep();
+       mMinTimestep = mpParam->getTimestep();
+       mMaxTimestep = mpParam->getTimestep();
 
        // init isosurf
-       this->mpIso->setIsolevel( this->mIsoValue );
+       mpIso->setIsolevel( mIsoValue );
 #if LBM_INCLUDE_TESTSOLVERS==1
        if(mUseTestdata) {
-               mpTest->setMaterialName( this->mpIso->getMaterialName() );
-               delete this->mpIso;
-               this->mpIso = mpTest;
-               if(mpTest->mDebugvalue1>0.0) { // 3d off
+               mpTest->setMaterialName( mpIso->getMaterialName() );
+               delete mpIso;
+               mpIso = mpTest;
+               if(mpTest->mFarfMode>0) { // 3d off
                        mpTest->setIsolevel(-100.0);
                } else {
-                       mpTest->setIsolevel( this->mIsoValue );
+                       mpTest->setIsolevel( mIsoValue );
                }
        }
-#endif // ELBEEM_PLUGIN!=1
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
        // approximate feature size with mesh resolution
        float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
        // smooth vars defined in solver_interface, set by simulation object
        // reset for invalid values...
-       if((this->mSmoothSurface<0.)||(this->mSmoothSurface>50.)) this->mSmoothSurface = 1.;
-       if((this->mSmoothNormals<0.)||(this->mSmoothNormals>50.)) this->mSmoothNormals = 1.;
-       this->mpIso->setSmoothSurface( this->mSmoothSurface * featureSize );
-       this->mpIso->setSmoothNormals( this->mSmoothNormals * featureSize );
+       if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
+       if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
+       mpIso->setSmoothSurface( mSmoothSurface * featureSize );
+       mpIso->setSmoothNormals( mSmoothNormals * featureSize );
 
        // init iso weight values mIsoWeightMethod
        int wcnt = 0;
@@ -835,34 +882,54 @@ bool LbmFsgrSolver::initializeSolverMemory()
                        }
        for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
 
-       LbmVec isostart = vec2L(this->mvGeoStart);
-       LbmVec isoend   = vec2L(this->mvGeoEnd);
+       LbmVec isostart = vec2L(mvGeoStart);
+       LbmVec isoend   = vec2L(mvGeoEnd);
        int twodOff = 0; // 2d slices
        if(LBMDIM==2) {
                LbmFloat sn,se;
-               sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5;
-               se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5;
+               sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
+               se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
                isostart[2] = sn;
                isoend[2] = se;
                twodOff = 2;
        }
-       int isosx = this->mSizex+2;
-       int isosy = this->mSizey+2;
-       int isosz = this->mSizez+2+twodOff;
+       int isosx = mSizex+2;
+       int isosy = mSizey+2;
+       int isosz = mSizez+2+twodOff;
 
        // MPT
 #if LBM_INCLUDE_TESTSOLVERS==1
-       if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
-               int xfac = 2;
-               isosx *= xfac;
-               isoend[0] = isostart[0] + (isoend[0]-isostart[0])*(LbmFloat)xfac;
+       //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
+       if( (mMpNum>0) && (mMpIndex==0) ) {
+               //? mpindex==0
+               // restore original value for node0
+               isosx       = mOrgSizeX + 2;
+               isostart[0] = mOrgStartX;
+               isoend[0]   = mOrgEndX;
        }
+       errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
+                       << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
+                       << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
+       errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
+                       << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
+                       << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 
-       //errMsg(" SETISO ", " "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5)<<" "<<(LbmFloat)(this->mSizex+1.0)<<" " );
+       errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
+       errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
+                       " grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
        mpIso->setStart( vec2G(isostart) );
        mpIso->setEnd(   vec2G(isoend) );
        LbmVec isodist = isoend-isostart;
+
+       int isosubs = mIsoSubdivs;
+       if(mFarFieldSize>1.) {
+               errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
+               isosubs = 1;
+               mpIso->setUseFulledgeArrays(true);
+       }
+       mpIso->setSubdivs(isosubs);
+
        mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
 
        // reset iso field
@@ -876,7 +943,7 @@ bool LbmFsgrSolver::initializeSolverMemory()
        for(int lev=0; lev<=mMaxRefine; lev++) {
                FSGR_FORIJK_BOUNDS(lev) {
                        RFLAG(lev,i,j,k,0) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
-                       if(!this->mAllfluid) {
+                       if(!mAllfluid) {
                                initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
                        } else {
                                initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); 
@@ -884,38 +951,43 @@ bool LbmFsgrSolver::initializeSolverMemory()
                }
        }
 
+
        if(LBMDIM==2) {
-               if(this->mOutputSurfacePreview) {
+               if(mOutputSurfacePreview) {
                        errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
-                       this->mOutputSurfacePreview = 0; }
+                       mOutputSurfacePreview = 0; }
+       }
+       if((glob_mpactive) && (glob_mpindex>0)) {
+               mOutputSurfacePreview = 0;
        }
+
 #if LBM_USE_GUI==1
-       if(this->mOutputSurfacePreview) {
+       if(mOutputSurfacePreview) {
                errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
-               this->mOutputSurfacePreview = 0; }
+               mOutputSurfacePreview = 0; }
 #endif // LBM_USE_GUI==1
-       if(this->mOutputSurfacePreview) {
+       if(mOutputSurfacePreview) {
                // same as normal one, but use reduced size
-               mpPreviewSurface = new IsoSurface( this->mIsoValue );
+               mpPreviewSurface = new IsoSurface( mIsoValue );
                mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
-               mpPreviewSurface->setIsolevel( this->mIsoValue );
+               mpPreviewSurface->setIsolevel( mIsoValue );
                // usually dont display for rendering
                mpPreviewSurface->setVisible( false );
 
                mpPreviewSurface->setStart( vec2G(isostart) );
                mpPreviewSurface->setEnd(   vec2G(isoend) );
                LbmVec pisodist = isoend-isostart;
-               LbmFloat pfac = this->mPreviewFactor;
-               mpPreviewSurface->initializeIsosurface( (int)(pfac*this->mSizex)+2, (int)(pfac*this->mSizey)+2, (int)(pfac*this->mSizez)+2, vec2G(pisodist) );
-               //mpPreviewSurface->setName( this->getName() + "preview" );
+               LbmFloat pfac = mPreviewFactor;
+               mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
+               //mpPreviewSurface->setName( getName() + "preview" );
                mpPreviewSurface->setName( "preview" );
        
-               debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*this->mSizex)<<","<<(pfac*this->mSizey)<<","<<(pfac*this->mSizez)<<" enabled",10);
+               debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
        }
 
        // init defaults
        mAvgNumUsedCells = 0;
-       this->mFixMass= 0.0;
+       mFixMass= 0.0;
        return true;
 }
 
@@ -933,23 +1005,23 @@ bool LbmFsgrSolver::initializeSolverGrids() {
 
        CellFlagType domainBoundType = CFInvalid;
        // TODO use normal object types instad...
-       if(this->mDomainBound.find(string("free")) != string::npos) {
+       if(mDomainBound.find(string("free")) != string::npos) {
                domainBoundType = CFBnd | CFBndFreeslip;
-       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<this->mDomainBound,10);
-       } else if(this->mDomainBound.find(string("part")) != string::npos) {
+       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
+       } else if(mDomainBound.find(string("part")) != string::npos) {
                domainBoundType = CFBnd | CFBndPartslip; // part slip type
-       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<this->mDomainPartSlipValue<<"), value:"<<this->mDomainBound,10);
+       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
        } else { 
                domainBoundType = CFBnd | CFBndNoslip;
-       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<this->mDomainBound,10);
+       debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
        }
 
        // use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
-       int domainobj = (int)(this->mpGiObjects->size());
+       int domainobj = (int)(mpGiObjects->size());
        domainBoundType |= (domainobj<<24);
        //for(int i=0; i<(int)(domainobj+0); i++) {
-               //errMsg("GEOIN","i"<<i<<" "<<(*this->mpGiObjects)[i]->getName());
-               //if((*this->mpGiObjects)[i] == this->mpIso) { //check...
+               //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
+               //if((*mpGiObjects)[i] == mpIso) { //check...
                //} 
        //}
        //errMsg("GEOIN"," dm "<<(domainBoundType>>24));
@@ -1050,19 +1122,19 @@ bool LbmFsgrSolver::initializeSolverGrids() {
        mInitialMass = 0.0;
        int inmCellCnt = 0;
        FSGR_FORIJK1(mMaxRefine) {
-               if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid) ) {
+               if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
                        LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); 
                        FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
                        mInitialMass += fluidRho;
                        inmCellCnt ++;
-               } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
+               } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
                        mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
                        inmCellCnt ++;
                }
        }
        mCurrentVolume = mCurrentMass = mInitialMass;
 
-       ParamVec cspv = this->mpParam->calculateCellSize();
+       ParamVec cspv = mpParam->calculateCellSize();
        if(LBMDIM==2) cspv[2] = 1.0;
        inmCellCnt = 1;
        double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
@@ -1071,7 +1143,7 @@ bool LbmFsgrSolver::initializeSolverGrids() {
 
        //mStartSymm = false;
 #if ELBEEM_PLUGIN!=1
-       if((LBMDIM==2)&&(this->mSizex<200)) {
+       if((LBMDIM==2)&&(mSizex<200)) {
                if(!checkSymmetry("init")) {
                        errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
                } else {
@@ -1094,18 +1166,18 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
                adaptGrid(lev);
                coarseRestrictFromFine(lev);
        }
-       this->markedClearList();
+       markedClearList();
        myTime_t fsgrtend = getTime(); 
-       if(!this->mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
+       if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
        mNumFsgrChanges = 0;
 
-       for(int l=0; l<this->cDirNum; l++) { 
+       for(int l=0; l<cDirNum; l++) { 
                LbmFloat area = 0.5 * 0.5 *0.5;
                if(LBMDIM==2) area = 0.5 * 0.5;
 
-               if(this->dfVecX[l]!=0) area *= 0.5;
-               if(this->dfVecY[l]!=0) area *= 0.5;
-               if(this->dfVecZ[l]!=0) area *= 0.5;
+               if(dfVecX[l]!=0) area *= 0.5;
+               if(dfVecY[l]!=0) area *= 0.5;
+               if(dfVecZ[l]!=0) area *= 0.5;
                mFsgrCellArea[l] = area;
        } // l
 
@@ -1123,6 +1195,8 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
        stepMain();
 
        // prepare once...
+       mpIso->setParticles(mpParticles, mPartDropMassSub);
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
        prepareVisualization();
        // copy again for stats counting
        for(int lev=0; lev<=mMaxRefine; lev++) {
@@ -1132,9 +1206,9 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
 
 
        // now really done...
-  debugOut("LbmFsgrSolver::initialize : SurfaceGen: SmsOrg("<<this->mSmoothSurface<<","<<this->mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<this->mpIso->getSmoothSurface()<<","<<this->mpIso->getSmoothNormals()<<") ",10);
-  debugOut("LbmFsgrSolver::initialize : Init done ... ",10);
-       this->mInitDone = 1;
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
+  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
+       mInitDone = 1;
 
 #if LBM_INCLUDE_TESTSOLVERS==1
        initCpdata();
@@ -1201,8 +1275,8 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
                                        const LbmVec oldov=objvel; /*DEBUG*/ \
                                        /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
                                        const LbmFloat partv = mObjectPartslips[OId]; \
-                                       /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, this->dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
-                                       /* m[l] = (RAC(ccel, this->dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
+                                       /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
+                                       /* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
                                        objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
                                }
 
@@ -1221,7 +1295,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        const int otherSet = mLevel[level].setOther;
        LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
        // for debugging - check targetTime check during DEFAULT STREAM
-       LbmFloat targetTime = mSimulationTime + this->mpParam->getTimestep();
+       LbmFloat targetTime = mSimulationTime + mpParam->getTimestep();
        if(mLastSimTime == targetTime) {
                debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
                return;
@@ -1236,16 +1310,15 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
        CellFlagType ntype = CFInvalid;
        // WARNING - copied from geo init!
-       int numobjs = (int)(this->mpGiObjects->size());
+       int numobjs = (int)(mpGiObjects->size());
        ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
        // 2d display as rectangles
        ntlVec3Gfx iniPos(0.0);
        if(LBMDIM==2) {
                dvec[2] = 1.0; 
-               iniPos = (this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))-(dvec*0.0);
+               iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
        } else {
-               iniPos = (this->mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
-               //iniPos[2] = this->mvGeoStart[2] + dvec[2]*getForZMin1();
+               iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
        }
        
        if( (int)mObjectMassMovnd.size() < numobjs) {
@@ -1258,21 +1331,20 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
        int nbored;
        for(int OId=0; OId<numobjs; OId++) {
-               ntlGeometryObject *obj = (*this->mpGiObjects)[OId];
+               ntlGeometryObject *obj = (*mpGiObjects)[OId];
                bool skip = false;
-               if(obj->getGeoInitId() != this->mLbmInitId) skip=true;
+               if(obj->getGeoInitId() != mLbmInitId) skip=true;
                if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
                if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
-               debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<this->mLbmInitId, 10);
                if(skip) continue;
+               debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);
 
                if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || 
                                (obj->getGeoInitType()&FGI_FLUID) && staticInit ) {
 
                        otype = ntype = CFInvalid;
                        switch(obj->getGeoInitType()) {
-                                       /*
-                               case FGI_BNDPART: 
+                               /* case FGI_BNDPART: // old, use noslip for moving part/free objs
                                case FGI_BNDFREE: 
                                        if(!staticInit) {
                                                errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
@@ -1311,8 +1383,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                        if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;
 
                        // copied from  recalculateObjectSpeeds
-                       mObjectSpeeds[OId] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
-                       debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
+                       mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
+                       debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
 
                        //vector<ntlVec3Gfx> tNormals;
                        vector<ntlVec3Gfx> *pNormals = NULL;
@@ -1324,7 +1396,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                // do two full update
                                // TODO tNormals handling!?
                                mMOIVerticesOld.clear();
-                               obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
+                               obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
                                monTrafo += mMOIVerticesOld.size();
                                obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
                                monTrafo += mMOIVertices.size();
@@ -1352,7 +1424,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                        // FIXME?
                                        if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
                                        // get old type - CHECK FIXME , timestep could have changed - cause trouble?
-                                       ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - this->mpParam->getTimestep(),sourceTime);
+                                       ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
                                        if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
                                }
                                if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
@@ -1360,6 +1432,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                LbmFloat massCheck = 0.;
                                int massReinits=0;                              
                                bool fillCells = (mObjectMassMovnd[OId]<=-1.);
+                               LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
+                               
 
                                // first pass - set new obs. cells
                                if(active) {
@@ -1397,30 +1471,19 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                        // compute fluid acceleration
                                                        FORDF1 {
                                                                if(rflagnb[l]&(CFFluid|CFInter)) { 
-                                                                       const LbmFloat ux = this->dfDvecX[l]*objvel[0];
-                                                                       const LbmFloat uy = this->dfDvecY[l]*objvel[1];
-                                                                       const LbmFloat uz = this->dfDvecZ[l]*objvel[2];
-
-                                                                       /*LbmFloat *rhodstCell = RACPNT(level, i+dfVecX[l],j+dfVecY[l],k+dfVecZ[l],workSet);
-                                                                       LbmFloat rhodst = RAC(rhodstCell,dC);
-                                                                       for(int ll=1; ll<19; ll++) { rhodst += RAC(rhodstCell,ll); }
-                                                                       rhodst = 1.; // DEBUG TEST! */
-
-                                                                       LbmFloat factor = 2. * this->dfLength[l] * 3.0 * (ux+uy+uz); // 
-                                                                               /* FSTEST 
-                                                                       */
+                                                                       const LbmFloat ux = dfDvecX[l]*objvel[0];
+                                                                       const LbmFloat uy = dfDvecY[l]*objvel[1];
+                                                                       const LbmFloat uz = dfDvecZ[l]*objvel[2];
+
+                                                                       LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // 
                                                                        if(ntype&(CFBndFreeslip|CFBndPartslip)) {
                                                                                // missing, diag mass checks...
-                                                                               //if(l<=LBMDIM*2) massCheck += factor;
-
                                                                                //if(l<=LBMDIM*2) factor *= 1.4142;
-                                                                               //if(l>LBMDIM*2) factor = 0.;
-                                                                               //else 
-                                                                               //factor = 0.; // TODO, optimize
                                                                                factor *= 2.0; // TODO, optimize
                                                                        } else {
                                                                                factor *= 1.2; // TODO, optimize
                                                                        }
+                                                                       factor *= impactCorrFactor; // use manual tweaking channel
                                                                        RAC(dstCell,l) = factor;
                                                                        massCheck += factor;
                                                                } else {
@@ -1461,28 +1524,19 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                        } 
                                                        CellFlagType settype = CFInvalid;
                                                        if(nbored&CFFluid) {
-                                                       //if(nbored&(CFFluid|CFInter)) {
                                                                settype = CFInter|CFNoInterpolSrc; 
                                                                rhomass = 1.5;
                                                                if(!fillCells) rhomass = 0.;
-                                                               //settype = CFFluid|CFNoInterpolSrc; rhomass=1.; // test
-                                                               //rhomass = 1.01; // DEBUGT
 
                                                                OBJVEL_CALC;
                                                                if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
-                                                               //settype=CFFluid; // rhomass = 1.; objvel = LbmVec(0.); // DEBUG TEST
 
                                                                // new interpolate values
                                                                LbmFloat avgrho = 0.0;
                                                                LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
                                                                interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
-                                                               //objvel = LbmVec(avgux,avguy,avguz); 
-                                                               //avgrho=1.;
                                                                initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
                                                                //errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
-
-                                                               // old - fixed init
-                                                               //initVelocityCell(level, i,j,k, settype, 1., rhomass, objvel );
                                                                massCheck += rhomass;
                                                        } else {
                                                                settype = CFEmpty; rhomass = 0.0;
@@ -1495,7 +1549,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                }  // wasactive
 
                                // only compute mass transfer when init is done
-                               if(this->mInitDone) {
+                               if(mInitDone) {
                                        errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
                                                " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass ); 
                                        mObjectMassMovnd[OId] += massCheck;
@@ -1512,8 +1566,6 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
                                                if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
                                                        //changeFlag(level, i,j,k, workSet, setflflag);
-                                                       //QCELL(level, i,j,k, workSet,dMass) = 1.; QCELL(level, i,j,k, workSet,dFfrac) = 1.;
-                                                       //initVelocityCell(level, i,j,k, setflflag, 1., 1., speed );
                                                        initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
                                                } 
                                                //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
@@ -1628,7 +1680,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        
 #undef POS2GRID_CHECK
        myTime_t monend = getTime();
-       debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
+       if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
        mLastSimTime = targetTime;
 }
 
@@ -1651,20 +1703,20 @@ bool LbmFsgrSolver::initGeometryFlags() {
        myTime_t geotimestart = getTime(); 
        ntlGeometryObject *pObj;
        ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
-       debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< this->mLbmInitId <<") v"<<dvec,3);
+       debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
        // WARNING - copied to movobj init!
 
        /* init object velocities, this has always to be called for init */
-       this->initGeoTree();
-       if(this->mAllfluid) { 
-               this->freeGeoTree();
+       initGeoTree();
+       if(mAllfluid) { 
+               freeGeoTree();
                return true; }
 
        // make sure moving obstacles are inited correctly
        // preinit moving obj points
-       int numobjs = (int)(this->mpGiObjects->size());
+       int numobjs = (int)(mpGiObjects->size());
        for(int o=0; o<numobjs; o++) {
-               ntlGeometryObject *obj = (*this->mpGiObjects)[o];
+               ntlGeometryObject *obj = (*mpGiObjects)[o];
                //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
                if(
                                ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
@@ -1677,20 +1729,20 @@ bool LbmFsgrSolver::initGeometryFlags() {
        }
 
        // max speed init
-       ntlVec3Gfx maxMovobjVelRw = this->getGeoMaxMovementVelocity( mSimulationTime, this->mpParam->getTimestep() );
-       ntlVec3Gfx maxMovobjVel = vec2G( this->mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
-       this->mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
-       LbmFloat allowMax = this->mpParam->getTadapMaxSpeed();  // maximum allowed velocity
+       ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
+       ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
+       mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
+       LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
        debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
-       if(this->mpParam->getSimulationMaxSpeed() > allowMax) {
+       if(mpParam->getSimulationMaxSpeed() > allowMax) {
                // similar to adaptTimestep();
-               LbmFloat nextmax = this->mpParam->getSimulationMaxSpeed();
-               LbmFloat newdt = this->mpParam->getTimestep() * (allowMax / nextmax); // newtr
-               debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< this->mpParam->getTimestep() <<" ",5);
-               this->mpParam->setDesiredTimestep( newdt );
-               this->mpParam->calculateAllMissingValues( mSimulationTime, this->mSilent );
-               maxMovobjVel = vec2G( this->mpParam->calculateLattVelocityFromRw( vec2P(this->getGeoMaxMovementVelocity(
-                                     mSimulationTime, this->mpParam->getTimestep() )) ));
+               LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
+               LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
+               debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
+               mpParam->setDesiredTimestep( newdt );
+               mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
+               maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
+                                     mSimulationTime, mpParam->getTimestep() )) ));
                debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
        }
        recalculateObjectSpeeds();
@@ -1710,24 +1762,25 @@ bool LbmFsgrSolver::initGeometryFlags() {
        // 2d display as rectangles
        if(LBMDIM==2) {
                dvec[2] = 0.0; 
-               iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))+(dvec*0.5);
+               iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
                //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
        } else {
-               iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
+               iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
        }
 
 
        // first init boundary conditions
        // invalid cells are set to empty afterwards
+       // TODO use floop macros!?
        for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
                for(int j=1;j<mLevel[level].lSizey-1;j++) {
                        for(int i=1;i<mLevel[level].lSizex-1;i++) {
                                ntype = CFInvalid;
                                
                                GETPOS(i,j,k);
-                               const bool inside = this->geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
+                               const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
                                if(inside) {
-                                       pObj = (*this->mpGiObjects)[OId];
+                                       pObj = (*mpGiObjects)[OId];
                                        switch( pObj->getGeoInitType() ){
                                        case FGI_MBNDINFLOW:  
                                          if(! pObj->getIsAnimated() ) {
@@ -1795,7 +1848,7 @@ bool LbmFsgrSolver::initGeometryFlags() {
                                ntype = CFInvalid;
                                int inits = 0;
                                GETPOS(i,j,k);
-                               const bool inside = this->geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
+                               const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
                                if(inside) {
                                        ntype = CFFluid;
                                }
@@ -1838,13 +1891,15 @@ bool LbmFsgrSolver::initGeometryFlags() {
                }
        }
 
-       this->freeGeoTree();
+       freeGeoTree();
        myTime_t geotimeend = getTime(); 
        debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 ); 
        //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
        return true;
 }
 
+#undef GETPOS
+
 /*****************************************************************************/
 /* init part for all freesurface testcases */
 void LbmFsgrSolver::initFreeSurfaces() {
@@ -1853,10 +1908,10 @@ void LbmFsgrSolver::initFreeSurfaces() {
 
        // set interface cells 
        FSGR_FORIJK1(mMaxRefine) {
-               if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid )) {
+               if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
                        FORDF1 {
-                               int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
-                               if( TESTFLAG( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr), CFEmpty ) ) {
+                               int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
+                               if( ( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr)& CFEmpty ) ) {
                                        LbmFloat arho=0., aux=0., auy=0., auz=0.;
                                        interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
                                        //errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
@@ -1870,13 +1925,13 @@ void LbmFsgrSolver::initFreeSurfaces() {
        // remove invalid interface cells 
        FSGR_FORIJK1(mMaxRefine) {
                // remove invalid interface cells 
-               if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
+               if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
                        int delit = 0;
                        int NBs = 0; // neighbor flags or'ed 
                        int noEmptyNB = 1;
 
                        FORDF1 {
-                               if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l ), CFEmpty ) ) {
+                               if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
                                        noEmptyNB = 0;
                                }
                                NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
@@ -1941,12 +1996,12 @@ void LbmFsgrSolver::initFreeSurfaces() {
 #endif // COMPRESSGRIDS==0
                for(int j=1;j<mLevel[lev].lSizey-1;++j) {
                for(int i=1;i<mLevel[lev].lSizex-1;++i) {
-                       if( TESTFLAG( RFLAG(lev, i,j,k, mLevel[lev].setCurr), CFInter) ) {
+                       if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
                                LbmFloat mass = 0.0;
                                //LbmFloat nbdiv;
                                //FORDF0 {
-                               for(int l=0;(l<this->cDirNum); l++) { 
-                                       int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+                               for(int l=0;(l<cDirNum); l++) { 
+                                       int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
                                        if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
                                                mass += 1.0;
                                        }
@@ -1957,7 +2012,7 @@ void LbmFsgrSolver::initFreeSurfaces() {
                                }
 
                                //errMsg(" I ", PRINT_IJK<<" m"<<mass );
-                               QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)this->cDirNum) );
+                               QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
                                QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
                        }
                }}}
@@ -2186,7 +2241,7 @@ bool LbmFsgrSolver::checkSymmetry(string idstring)
        LbmFloat maxdiv =0.0;
        if(erro) {
                errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
-               //if(LBMDIM==2) this->mPanic = true; 
+               //if(LBMDIM==2) mPanic = true; 
                //return false;
        } else {
                errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );