- New options for mesh voxelization: shell only (also
authorNils Thuerey <nils@thuerey.de>
Thu, 11 May 2006 08:09:02 +0000 (08:09 +0000)
committerNils Thuerey <nils@thuerey.de>
Thu, 11 May 2006 08:09:02 +0000 (08:09 +0000)
  works for non closed objects), volume ("normal"/old way of
  doing it), and a combination of both:
  http://www10.informatik.uni-erlangen.de/~sinithue/blender/voltcomp_sm.jpg
- Finally included bjornmose MSVC6 fixes
- Added support for animated meshes, e.g. meshes with
  parented skeletons. Is enabled for obstacles with a new button.
  A simple example with Bassam's mancandy can be found here:
  http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid2_mancandy.mpg
  http://www10.informatik.uni-erlangen.de/~sinithue/blender/fluid2_mancandy.blend
  (Warning - keep meshes as simple as possible, e.g. turn off subsurf
  for baking. Export probably shoulb be further optimized.)
- Changed handling of no/free/part slip obstacles, see:
  http://www10.informatik.uni-erlangen.de/~sinithue/blender/bndtcomp_sm.jpg
- Removed surface particle option for upcoming release,
  needs more testing & tweaking
- Added tracer particles instead (swimming along in the fluid)
- Updated wiki (description of IPOs still missing).

41 files changed:
intern/elbeem/extern/LBM_fluidsim.h
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/ntl_blenderdumper.cpp
intern/elbeem/intern/ntl_blenderdumper.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_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 [new file with mode: 0644]
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/makesdna/DNA_object_fluidsim.h
source/blender/src/buttons_object.c
source/blender/src/fluidsim.c

index c93d129..a6ba16a 100644 (file)
@@ -52,7 +52,7 @@ void fluidsimSettingsFree(struct FluidsimSettings* sb);
 void fluidsimBake(struct Object* ob);
 
 /* read & write bobj / bobj.gz files (e.g. for fluid sim surface meshes) */
-void writeBobjgz(char *filename, struct Object *ob);
+void writeBobjgz(char *filename, struct Object *ob, int useGlobalCoords, int append, float time);
 struct Mesh* readBobjgz(char *filename, struct Mesh *orgmesh, float* bbstart, float *bbsize);
 
 /* create derived mesh for fluid sim objects */
index f0d154c..8eb8a54 100644 (file)
 #define ELBEEM_API_H
 
 
-/*! blender types for mesh->type */
-//#define OB_FLUIDSIM_DOMAIN      2
-#define OB_FLUIDSIM_FLUID       4
-#define OB_FLUIDSIM_OBSTACLE    8
-#define OB_FLUIDSIM_INFLOW      16
-#define OB_FLUIDSIM_OUTFLOW     32
-
-#define FLUIDSIM_OBSTACLE_NOSLIP     1
-#define FLUIDSIM_OBSTACLE_PARTSLIP   2
-#define FLUIDSIM_OBSTACLE_FREESLIP   3
+// simulation run callback function type (elbeemSimulationSettings->runsimCallback)
+// best use with FLUIDSIM_CBxxx defines below.
+// >parameters
+// return values: 0=continue, 1=stop, 2=abort
+// data pointer: user data pointer from elbeemSimulationSettings->runsimUserData
+// status integer: 1=running simulation, 2=new frame saved
+// frame integer: if status is 1, contains current frame number
+typedef int (*elbeemRunSimulationCallback)(void *data, int status, int frame);
+#define FLUIDSIM_CBRET_CONTINUE    0
+#define FLUIDSIM_CBRET_STOP        1
+#define FLUIDSIM_CBRET_ABORT       2 
+#define FLUIDSIM_CBSTATUS_STEP     1 
+#define FLUIDSIM_CBSTATUS_NEWFRAME 2 
 
 
 // global settings for the simulation
 typedef struct elbeemSimulationSettings {
   /* version number */
   short version;
+       /* id number of simulation domain, needed if more than a
+        * single domain should be simulated */
+       short domainId;
 
        /* geometrical extent */
        float geoStart[3], geoSize[3];
@@ -49,8 +55,10 @@ typedef struct elbeemSimulationSettings {
   float gstar;
   /* activate refinement? */
   short maxRefine;
-  /* amount of particles to generate (0=off) */
+  /* probability for surface particle generation (0.0=off) */
   float generateParticles;
+  /* amount of tracer particles to generate (0=off) */
+  int numTracerParticles;
 
   /* store output path, and file prefix for baked fluid surface */
   char outputPath[160+80];
@@ -77,17 +85,43 @@ typedef struct elbeemSimulationSettings {
        /* development variables, testing for upcoming releases...*/
        float farFieldSize;
 
+       /* callback function to notify calling program of performed simulation steps
+        * or newly available frame data, if NULL it is ignored */
+       elbeemRunSimulationCallback runsimCallback;
+       /* pointer passed to runsimCallback for user data storage */
+       void* runsimUserData;
+
 } elbeemSimulationSettings;
 
 
+// defines for elbeemMesh->type below
+#define OB_FLUIDSIM_FLUID       4
+#define OB_FLUIDSIM_OBSTACLE    8
+#define OB_FLUIDSIM_INFLOW      16
+#define OB_FLUIDSIM_OUTFLOW     32
+
+// defines for elbeemMesh->obstacleType below
+#define FLUIDSIM_OBSTACLE_NOSLIP     1
+#define FLUIDSIM_OBSTACLE_PARTSLIP   2
+#define FLUIDSIM_OBSTACLE_FREESLIP   3
+
+#define OB_VOLUMEINIT_VOLUME 1
+#define OB_VOLUMEINIT_SHELL  2
+#define OB_VOLUMEINIT_BOTH   (OB_VOLUMEINIT_SHELL|OB_VOLUMEINIT_VOLUME)
+
 // a single mesh object
 typedef struct elbeemMesh {
   /* obstacle,fluid or inflow... */
   short type;
+       /* id of simulation domain it belongs to */
+       short parentDomainId;
 
        /* vertices */
   int numVertices;
-       float *vertices; // = float[][3];
+       float *vertices; // = float[n][3];
+       /* animated vertices */
+  int channelSizeVertices;
+       float *channelVertices; // = float[channelSizeVertices* (n*3+1) ];
 
        /* triangles */
        int   numTriangles;
@@ -112,6 +146,8 @@ typedef struct elbeemMesh {
        /* boundary types and settings */
        short obstacleType;
        float obstaclePartslip;
+       /* init volume, shell or both? use OB_VOLUMEINIT_xxx defines above */
+       short volumeInitType;
 
        /* name of the mesh, mostly for debugging */
        char *name;
@@ -123,11 +159,16 @@ typedef struct elbeemMesh {
 extern "C"  {
 #endif // __cplusplus
  
+
 // reset elbeemSimulationSettings struct with defaults
 void elbeemResetSettings(struct elbeemSimulationSettings*);
  
 // start fluidsim init (returns !=0 upon failure)
-int elbeemInit(struct elbeemSimulationSettings*);
+int elbeemInit(void);
+
+// start fluidsim init (returns !=0 upon failure)
+int elbeemAddDomain(struct elbeemSimulationSettings*);
+
 // get failure message during simulation or init
 // if an error occured (the string is copied into buffer,
 // max. length = 256 chars )
@@ -142,6 +183,9 @@ int elbeemAddMesh(struct elbeemMesh*);
 // do the actual simulation
 int elbeemSimulate(void);
 
+// continue a previously stopped simulation
+int elbeemContinueSimulation(void);
+
 
 // helper functions 
 
@@ -185,7 +229,5 @@ double elbeemEstimateMemreq(int res,
 #define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6))
 #define FGI_MBNDOUTFLOW        (1<<(FGI_FLAGSTART+ 7))
 
-#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
-
 
 #endif // ELBEEM_API_H
index abdd931..f970f6e 100644 (file)
@@ -119,7 +119,7 @@ int Attribute::getAsInt()
                errMsg("Attribute::getAsInt", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
                errMsg("Attribute::getAsInt", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
 #if ELBEEM_PLUGIN!=1
-               gElbeemState = -4; // parse error
+               setElbeemState( -4 ); // parse error
 #endif
                return 0;
        }
@@ -159,7 +159,7 @@ double Attribute::getAsFloat()
                errMsg("Attribute::getAsFloat", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
                errMsg("Attribute::getAsFloat", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
 #if ELBEEM_PLUGIN!=1
-               gElbeemState = -4; // parse error
+               setElbeemState( -4 ); // parse error
 #endif
                return 0.0;
        }
@@ -211,7 +211,7 @@ ntlVec3d Attribute::getAsVec3d()
                errMsg("Attribute::getAsVec3d", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
                errMsg("Attribute::getAsVec3d", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
 #if ELBEEM_PLUGIN!=1
-               gElbeemState = -4; // parse error
+               setElbeemState( -4 ); // parse error
 #endif
                return ntlVec3d(0.0);
        }
@@ -267,7 +267,7 @@ void Attribute::getAsMat4Gfx(ntlMat4Gfx *mat)
                errMsg("Attribute::getAsMat4Gfx", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
                errMsg("Attribute::getAsMat4Gfx", "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" );
 #if ELBEEM_PLUGIN!=1
-               gElbeemState = -4; // parse error
+               setElbeemState( -4 ); // parse error
 #endif
                *mat = ntlMat4Gfx(0.0);
                return;
@@ -349,6 +349,30 @@ AnimChannel<ntlVec3d> Attribute::getChannelVec3d() {
        return AnimChannel<ntlVec3d>(valVec,timeVec);
 }
 
+//! get channel as float vector set
+AnimChannel<ntlSetVec3f> 
+Attribute::getChannelSetVec3f() {
+       vector<double> timeVec;
+       ntlSetVec3f setv;
+       vector< ntlSetVec3f > valVec;
+       
+       if((!initChannel(3)) || (!mIsChannel)) {
+               timeVec.push_back( 0.0 );
+               setv.mVerts.push_back( vec2F( getAsVec3d() ) );
+               valVec.push_back( setv );
+       } else {
+               for(size_t i=0; i<mChannel.size(); i++) {
+                       mValue = mChannel[i];
+                       ntlVec3f val = vec2F( getAsVec3d() );
+                       timeVec.push_back( mTimes[i] );
+                       setv.mVerts.push_back( val );
+               }
+               valVec.push_back( setv );
+               valVec.push_back( setv );
+       }
+
+       return AnimChannel<ntlSetVec3f>(valVec,timeVec);
+}
 
 /******************************************************************************
  * check if there were unknown params
@@ -447,29 +471,50 @@ bool AttributeList::ignoreParameter(string name, string source) {
 }
                
 // read channels
-AnimChannel<double> AttributeList::readChannelFloat(string name) {
-       if(!exists(name)) { return AnimChannel<double>(0.0); } 
-       AnimChannel<double> ret = find(name)->getChannelFloat(); 
+AnimChannel<int> AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) {
+       if(!exists(name)) { return AnimChannel<int>(defaultValue); } 
+       AnimChannel<int> ret = find(name)->getChannelInt(); 
        find(name)->setUsed(true);
-       channelSimplifyd(ret);
+       channelSimplifyi(ret);
        return ret;
 }
-AnimChannel<int> AttributeList::readChannelInt(string name) {
-       if(!exists(name)) { return AnimChannel<int>(0); } 
-       AnimChannel<int> ret = find(name)->getChannelInt(); 
+AnimChannel<double> AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) {
+       if(!exists(name)) { return AnimChannel<double>(defaultValue); } 
+       AnimChannel<double> ret = find(name)->getChannelFloat(); 
        find(name)->setUsed(true);
-       channelSimplifyi(ret);
+       channelSimplifyd(ret);
        return ret;
 }
-AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name) {
-       if(!exists(name)) { return AnimChannel<ntlVec3d>(0.0); } 
+AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) {
+       if(!exists(name)) { return AnimChannel<ntlVec3d>(defaultValue); } 
        AnimChannel<ntlVec3d> ret = find(name)->getChannelVec3d(); 
        find(name)->setUsed(true);
        channelSimplifyVd(ret);
        return ret;
 }
-AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name) {
-       if(!exists(name)) { return AnimChannel<ntlVec3f>(0.0); } 
+AnimChannel<ntlSetVec3f> AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) {
+       if(!exists(name)) { return AnimChannel<ntlSetVec3f>(defaultValue); } 
+       AnimChannel<ntlSetVec3f> ret = find(name)->getChannelSetVec3f(); 
+       find(name)->setUsed(true);
+       //channelSimplifyVf(ret);
+       return ret;
+}
+AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) {
+       if(!exists(name)) { 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;
+}
+AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) {
+       if(!exists(name)) { return AnimChannel<ntlVec3f>(defaultValue); } 
 
        AnimChannel<ntlVec3d> convert = find(name)->getChannelVec3d(); 
        // convert to float
@@ -672,6 +717,7 @@ bool channelSimplifyVd (AnimChannel<ntlVec3d> &channel) {
        return channelSimplifyVecT<ntlVec3d>(channel);
 }
 
+//! debug function, prints channel as string
 template<class Scalar>
 string AnimChannel<Scalar>::printChannel() {
        std::ostringstream ostr;
@@ -684,5 +730,65 @@ 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() ); }
+}
+
+
+ntlSetVec3f::ntlSetVec3f(double v ) {
+       mVerts.clear();
+       mVerts.push_back( ntlVec3f(v) );
+}
+const ntlSetVec3f& 
+ntlSetVec3f::operator=(double v ) {
+       mVerts.clear();
+       mVerts.push_back( ntlVec3f(v) );
+       return *this;
+}
+
+std::ostream& operator<<( std::ostream& os, const ntlSetVec3f& vs ) {
+       os<< "{";
+       for(int j=0;j<(int)vs.mVerts.size();j++)  os<<vs.mVerts[j];
+       os<< "}";
+  return os;
+}
+
+ntlSetVec3f& 
+ntlSetVec3f::operator+=( double v )
+{
+       for(int j=0;j<(int)(mVerts.size()) ;j++) {
+               mVerts[j] += v;
+       }
+       return *this;
+}
+
+ntlSetVec3f& 
+ntlSetVec3f::operator+=( const ntlSetVec3f &v )
+{
+       for(int j=0;j<(int)MIN(mVerts.size(),v.mVerts.size()) ;j++) {
+               mVerts[j] += v.mVerts[j];
+       }
+       return *this;
+}
+
+ntlSetVec3f& 
+ntlSetVec3f::operator*=( double v )
+{
+       for(int j=0;j<(int)(mVerts.size()) ;j++) {
+               mVerts[j] *= v;
+       }
+       return *this;
+}
+
+ntlSetVec3f& 
+ntlSetVec3f::operator*=( const ntlSetVec3f &v )
+{
+       for(int j=0;j<(int)MIN(mVerts.size(),v.mVerts.size()) ;j++) {
+               mVerts[j] *= v.mVerts[j];
+       }
+       return *this;
+}
 
 
index 730e136..20e5341 100644 (file)
@@ -12,6 +12,9 @@
 
 #include "utilities.h"
 template<class T> class ntlMatrix4x4;
+class ntlSetVec3f;
+std::ostream& operator<<( std::ostream& os, const ntlSetVec3f& i );
+
 
 
 //! An animated attribute channel
@@ -21,15 +24,15 @@ class AnimChannel
        public:
                // default constructor
                AnimChannel() : 
-                       mValue(), mTimes() { mInited = false; }
+                       mValue(), mTimes() { mInited = false; debugPrintChannel(); }
 
                // null init constructor
                AnimChannel(Scalar null) : 
-                       mValue(1), mTimes(1) { mValue[0]=null; mTimes[0]=0.0; mInited = true; }
+                       mValue(1), mTimes(1) { mValue[0]=null; mTimes[0]=0.0; mInited = true; debugPrintChannel(); }
 
                // proper init
-               AnimChannel(vector<Scalar> v, vector<double> t) : 
-                       mValue(v), mTimes(t) { mInited = true; }
+               AnimChannel(vector<Scalar> &v, vector<double> &t) : 
+                       mValue(v), mTimes(t) { mInited = true; debugPrintChannel(); }
 
                // desctructor, nothing to do
                ~AnimChannel() { };
@@ -45,7 +48,14 @@ class AnimChannel
                                        // interpolate
                                        double d = mTimes[i+1]-mTimes[i];
                                        double f = (t-mTimes[i])/d;
-                                       return (Scalar)(mValue[i] * (1.0-f) + mValue[i+1] * f);
+                                       //return (Scalar)(mValue[i] * (1.0-f) + mValue[i+1] * f);
+                                       Scalar ret,tmp;
+                                       ret = mValue[i];
+                                       ret *= 1.-f;
+                                       tmp = mValue[i+1];
+                                       tmp *= f;
+                                       ret += tmp;
+                                       return ret;
                                }
                        }
                        // whats this...?
@@ -77,6 +87,8 @@ class AnimChannel
 
                //! debug function, prints channel as string
                string printChannel();
+               //! debug function, prints to stdout if DEBUG_CHANNELS flag is enabled, used in constructors
+               void debugPrintChannel();
                //! valid init?
                bool isInited() { return mInited; }
 
@@ -145,6 +157,8 @@ class Attribute
                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();
@@ -181,6 +195,22 @@ class 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;
+};
+
 //! The list of configuration attributes
 class AttributeList
 {
@@ -230,10 +260,13 @@ class AttributeList
                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);
-               AnimChannel<double>  readChannelFloat(string name);
-               AnimChannel<ntlVec3d> readChannelVec3d(string name);
-               AnimChannel<ntlVec3f> readChannelVec3f(string name);
+               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);
index daee5f8..bd96b9f 100644 (file)
@@ -36,7 +36,8 @@ ntlWorld *gpWorld = NULL;
 extern "C" 
 void elbeemResetSettings(elbeemSimulationSettings *set) {
        if(!set) return;
-  set->version = 2;
+  set->version = 3;
+  set->domainId = 0;
        for(int i=0 ; i<3; i++) set->geoStart[i] = 0.0;
        for(int i=0 ; i<3; i++) set->geoSize[i] = 1.0;
   set->resolutionxyz = 64;
@@ -52,7 +53,8 @@ void elbeemResetSettings(elbeemSimulationSettings *set) {
        set->noOfFrames = 10;
   set->gstar = 0.005;
   set->maxRefine = -1;
-  set->generateParticles = 1.0;
+  set->generateParticles = 0.0;
+  set->numTracerParticles = 0;
   strcpy(set->outputPath,"./elbeemdata_");
 
        set->channelSizeFrameTime=0;
@@ -68,6 +70,8 @@ void elbeemResetSettings(elbeemSimulationSettings *set) {
        set->surfaceSmoothing = 1.;
 
        set->farFieldSize = 0.;
+       set->runsimCallback = NULL;
+       set->runsimUserData = NULL;
 
        // init identity
        for(int i=0; i<16; i++) set->surfaceTrafo[i] = 0.0;
@@ -76,24 +80,35 @@ void elbeemResetSettings(elbeemSimulationSettings *set) {
 
 // start fluidsim init
 extern "C" 
-int elbeemInit(elbeemSimulationSettings *settings) {
-       gElbeemState = SIMWORLD_INVALID;
-       strcpy(gElbeemErrorString,"[none]");
+int elbeemInit() {
+       setElbeemState( SIMWORLD_INITIALIZING );
+       setElbeemErrorString("[none]");
 
        elbeemCheckDebugEnv();
        debMsgStd("performElbeemSimulation",DM_NOTIFY,"El'Beem Simulation Init Start as Plugin, debugLevel:"<<gDebugLevel<<" ...\n", 2);
        
        // create world object with initial settings
-       ntlBlenderDumper *elbeem = new ntlBlenderDumper(settings);
+       ntlBlenderDumper *elbeem = new ntlBlenderDumper(); 
        gpWorld = elbeem;
        return 0;
 }
 
+// start fluidsim init
+extern "C" 
+int elbeemAddDomain(elbeemSimulationSettings *settings) {
+       // has to be inited...
+       if((getElbeemState() == SIMWORLD_INVALID) && (!gpWorld)) { elbeemInit(); }
+       if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddDomain","Unable to init simulation world",SIMWORLD_INITERROR); }
+       // create domain with given settings
+       gpWorld->addDomain(settings);
+       return 0;
+}
+
 // error message access
 extern "C" 
 void elbeemGetErrorString(char *buffer) {
        if(!buffer) return;
-       strncpy(buffer,gElbeemErrorString,256);
+       strncpy(buffer,getElbeemErrorString(),256);
 }
 
 // reset elbeemMesh struct with zeroes
@@ -105,6 +120,7 @@ void elbeemResetMesh(elbeemMesh *mesh) {
        mesh->vertices = NULL;
        mesh->numTriangles = 0;
   mesh->triangles = NULL;
+
        mesh->channelSizeTranslation = 0;
        mesh->channelTranslation = NULL;
        mesh->channelSizeRotation = 0;
@@ -116,15 +132,24 @@ void elbeemResetMesh(elbeemMesh *mesh) {
        mesh->channelSizeInitialVel = 0;
        mesh->channelInitialVel = NULL;
        mesh->localInivelCoords = 0;
+
        mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP;
+       mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
        mesh->obstaclePartslip= 0.;
+
+       mesh->channelSizeVertices = 0;
+       mesh->channelVertices = NULL;
+
        mesh->name = "[unnamed]";
 }
 
+int globalMeshCounter = 1;
 // add mesh as fluidsim object
 extern "C" 
 int elbeemAddMesh(elbeemMesh *mesh) {
        int initType = -1;
+       if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddMesh","World and domain not initialized, call elbeemInit and elbeemAddDomain before...", SIMWORLD_INITERROR); }
+
        switch(mesh->type) {
                case OB_FLUIDSIM_OBSTACLE: 
                        if     (mesh->obstacleType==FLUIDSIM_OBSTACLE_PARTSLIP) initType = FGI_BNDPART; 
@@ -140,13 +165,22 @@ int elbeemAddMesh(elbeemMesh *mesh) {
        
        ntlGeometryObjModel *obj = new ntlGeometryObjModel( );
        gpWorld->getRenderGlobals()->getSimScene()->addGeoClass( obj );
-       obj->initModel(mesh->numVertices, mesh->vertices, mesh->numTriangles, mesh->triangles);
-       if(mesh->name) obj->setName(std::string(mesh->name));
-       else obj->setName(std::string("[unnamed]"));
-       obj->setGeoInitId(1);
+       obj->initModel(
+                       mesh->numVertices, mesh->vertices, mesh->numTriangles, mesh->triangles,
+                       mesh->channelSizeVertices, mesh->channelVertices );
+       if(mesh->name) obj->setName(string(mesh->name));
+       else {
+               char meshname[100];
+               snprintf(meshname,100,"mesh%04d",globalMeshCounter);
+               obj->setName(string(meshname));
+       }
+       globalMeshCounter++;
+       obj->setGeoInitId( mesh->parentDomainId+1 );
        obj->setGeoInitIntersect(true);
        obj->setGeoInitType(initType);
        obj->setGeoPartSlipValue(mesh->obstaclePartslip);
+       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]) );
        obj->initChannels(
                        mesh->channelSizeTranslation, mesh->channelTranslation, 
@@ -167,17 +201,19 @@ int elbeemSimulate(void) {
        if(!gpWorld) return 1;
 
        gpWorld->finishWorldInit();
-       
-       if(SIMWORLD_OK()) {
-               gElbeemState = SIMWORLD_INITED;
+       if( isSimworldOk() ) {
                myTime_t timestart = getTime();
                gpWorld->renderAnimation();
                myTime_t timeend = getTime();
-               debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
 
-               // ok, we're done...
-               delete gpWorld;
-               gpWorld = NULL;
+               if(getElbeemState() != SIMWORLD_STOP) {
+                       // ok, we're done...
+                       delete gpWorld;
+                       gpWorld = NULL;
+                       debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
+               } else {
+                       debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
+               }
                return 0;
        } 
 
@@ -186,6 +222,32 @@ int elbeemSimulate(void) {
 }
 
 
+// continue a previously stopped simulation
+extern "C" 
+int elbeemContinueSimulation(void) {
+
+       if(getElbeemState() != SIMWORLD_STOP) {
+               errMsg("elbeemContinueSimulation","No running simulation found! Aborting...");
+               if(gpWorld) delete gpWorld;
+               return 1;
+       }
+
+       myTime_t timestart = getTime();
+       gpWorld->renderAnimation();
+       myTime_t timeend = getTime();
+
+       if(getElbeemState() != SIMWORLD_STOP) {
+               // ok, we're done...
+               delete gpWorld;
+               gpWorld = NULL;
+               debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
+       } else {
+               debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
+       }
+       return 0;
+} 
+
+
 // global vector to flag values to remove
 vector<int> gKeepVal;
 
index f0d154c..8eb8a54 100644 (file)
 #define ELBEEM_API_H
 
 
-/*! blender types for mesh->type */
-//#define OB_FLUIDSIM_DOMAIN      2
-#define OB_FLUIDSIM_FLUID       4
-#define OB_FLUIDSIM_OBSTACLE    8
-#define OB_FLUIDSIM_INFLOW      16
-#define OB_FLUIDSIM_OUTFLOW     32
-
-#define FLUIDSIM_OBSTACLE_NOSLIP     1
-#define FLUIDSIM_OBSTACLE_PARTSLIP   2
-#define FLUIDSIM_OBSTACLE_FREESLIP   3
+// simulation run callback function type (elbeemSimulationSettings->runsimCallback)
+// best use with FLUIDSIM_CBxxx defines below.
+// >parameters
+// return values: 0=continue, 1=stop, 2=abort
+// data pointer: user data pointer from elbeemSimulationSettings->runsimUserData
+// status integer: 1=running simulation, 2=new frame saved
+// frame integer: if status is 1, contains current frame number
+typedef int (*elbeemRunSimulationCallback)(void *data, int status, int frame);
+#define FLUIDSIM_CBRET_CONTINUE    0
+#define FLUIDSIM_CBRET_STOP        1
+#define FLUIDSIM_CBRET_ABORT       2 
+#define FLUIDSIM_CBSTATUS_STEP     1 
+#define FLUIDSIM_CBSTATUS_NEWFRAME 2 
 
 
 // global settings for the simulation
 typedef struct elbeemSimulationSettings {
   /* version number */
   short version;
+       /* id number of simulation domain, needed if more than a
+        * single domain should be simulated */
+       short domainId;
 
        /* geometrical extent */
        float geoStart[3], geoSize[3];
@@ -49,8 +55,10 @@ typedef struct elbeemSimulationSettings {
   float gstar;
   /* activate refinement? */
   short maxRefine;
-  /* amount of particles to generate (0=off) */
+  /* probability for surface particle generation (0.0=off) */
   float generateParticles;
+  /* amount of tracer particles to generate (0=off) */
+  int numTracerParticles;
 
   /* store output path, and file prefix for baked fluid surface */
   char outputPath[160+80];
@@ -77,17 +85,43 @@ typedef struct elbeemSimulationSettings {
        /* development variables, testing for upcoming releases...*/
        float farFieldSize;
 
+       /* callback function to notify calling program of performed simulation steps
+        * or newly available frame data, if NULL it is ignored */
+       elbeemRunSimulationCallback runsimCallback;
+       /* pointer passed to runsimCallback for user data storage */
+       void* runsimUserData;
+
 } elbeemSimulationSettings;
 
 
+// defines for elbeemMesh->type below
+#define OB_FLUIDSIM_FLUID       4
+#define OB_FLUIDSIM_OBSTACLE    8
+#define OB_FLUIDSIM_INFLOW      16
+#define OB_FLUIDSIM_OUTFLOW     32
+
+// defines for elbeemMesh->obstacleType below
+#define FLUIDSIM_OBSTACLE_NOSLIP     1
+#define FLUIDSIM_OBSTACLE_PARTSLIP   2
+#define FLUIDSIM_OBSTACLE_FREESLIP   3
+
+#define OB_VOLUMEINIT_VOLUME 1
+#define OB_VOLUMEINIT_SHELL  2
+#define OB_VOLUMEINIT_BOTH   (OB_VOLUMEINIT_SHELL|OB_VOLUMEINIT_VOLUME)
+
 // a single mesh object
 typedef struct elbeemMesh {
   /* obstacle,fluid or inflow... */
   short type;
+       /* id of simulation domain it belongs to */
+       short parentDomainId;
 
        /* vertices */
   int numVertices;
-       float *vertices; // = float[][3];
+       float *vertices; // = float[n][3];
+       /* animated vertices */
+  int channelSizeVertices;
+       float *channelVertices; // = float[channelSizeVertices* (n*3+1) ];
 
        /* triangles */
        int   numTriangles;
@@ -112,6 +146,8 @@ typedef struct elbeemMesh {
        /* boundary types and settings */
        short obstacleType;
        float obstaclePartslip;
+       /* init volume, shell or both? use OB_VOLUMEINIT_xxx defines above */
+       short volumeInitType;
 
        /* name of the mesh, mostly for debugging */
        char *name;
@@ -123,11 +159,16 @@ typedef struct elbeemMesh {
 extern "C"  {
 #endif // __cplusplus
  
+
 // reset elbeemSimulationSettings struct with defaults
 void elbeemResetSettings(struct elbeemSimulationSettings*);
  
 // start fluidsim init (returns !=0 upon failure)
-int elbeemInit(struct elbeemSimulationSettings*);
+int elbeemInit(void);
+
+// start fluidsim init (returns !=0 upon failure)
+int elbeemAddDomain(struct elbeemSimulationSettings*);
+
 // get failure message during simulation or init
 // if an error occured (the string is copied into buffer,
 // max. length = 256 chars )
@@ -142,6 +183,9 @@ int elbeemAddMesh(struct elbeemMesh*);
 // do the actual simulation
 int elbeemSimulate(void);
 
+// continue a previously stopped simulation
+int elbeemContinueSimulation(void);
+
 
 // helper functions 
 
@@ -185,7 +229,5 @@ double elbeemEstimateMemreq(int res,
 #define FGI_MBNDINFLOW (1<<(FGI_FLAGSTART+ 6))
 #define FGI_MBNDOUTFLOW        (1<<(FGI_FLAGSTART+ 7))
 
-#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
-
 
 #endif // ELBEEM_API_H
index 07bc6b7..373a6f7 100644 (file)
@@ -303,7 +303,7 @@ void IsoSurface::triangulate( void )
 /******************************************************************************
  * Get triangles for rendering
  *****************************************************************************/
-void IsoSurface::getTriangles( vector<ntlTriangle> *triangles, 
+void IsoSurface::getTriangles(double t, vector<ntlTriangle> *triangles, 
                                                                                                         vector<ntlVec3Gfx> *vertices, 
                                                                                                         vector<ntlVec3Gfx> *normals, int objectId )
 {
@@ -311,6 +311,7 @@ void IsoSurface::getTriangles( vector<ntlTriangle> *triangles,
                debugOut("IsoSurface::getTriangles warning: Not initialized! ", 10);
                return;
        }
+       t = 0.;
        //return; // DEBUG
 
   /* triangulate field */
@@ -355,12 +356,6 @@ void IsoSurface::getTriangles( vector<ntlTriangle> *triangles,
                if(getVisible()){ flag |= TRI_GEOMETRY; }
                if(getCastShadows() ) { 
                        flag |= TRI_CASTSHADOWS; } 
-               if( (getMaterial()->getMirror()>0.0) ||  
-                               (getMaterial()->getTransparence()>0.0) ||  
-                               (getMaterial()->getFresnel()>0.0) ) { 
-                       flag |= TRI_MAKECAUSTICS; } 
-               else { 
-                       flag |= TRI_NOCAUSTICS; } 
 
                /* init geo init id */
                int geoiId = getGeoInitId(); 
index 5e5115e..d1ff152 100644 (file)
@@ -123,7 +123,7 @@ class IsoSurface :
                inline float getSmoothNormals() { return mSmoothNormals; }
 
                // geometry object functions 
-               virtual void getTriangles( vector<ntlTriangle> *triangles, 
+               virtual void getTriangles(double t, vector<ntlTriangle> *triangles, 
                                vector<ntlVec3Gfx> *vertices, 
                                vector<ntlVec3Gfx> *normals, int objectId );
 
index ecd1967..4468dda 100644 (file)
 /******************************************************************************
  * Constructor
  *****************************************************************************/
-ntlBlenderDumper::ntlBlenderDumper(string filename, bool commandlineMode) :
-       ntlWorld(filename,commandlineMode),
-       mpTrafo(NULL)
+ntlBlenderDumper::ntlBlenderDumper() : ntlWorld()
 {
-  ntlRenderGlobals *glob = mpGlob;
-       AttributeList *pAttrs = glob->getBlenderAttributes();
-       mpTrafo = new ntlMat4Gfx(0.0);
-       mpTrafo->initId();
-       pAttrs->readMat4Gfx("transform" , (*mpTrafo), "ntlBlenderDumper","mpTrafo", false, mpTrafo ); 
+       // same as normal constructor here
 }
-ntlBlenderDumper::ntlBlenderDumper(elbeemSimulationSettings *settings) :
-       ntlWorld(settings), mpTrafo(NULL)
+ntlBlenderDumper::ntlBlenderDumper(string filename, bool commandlineMode) :
+       ntlWorld(filename,commandlineMode)
 {
-       // same as normal constructor here
-       mpTrafo = new ntlMat4Gfx(0.0);
-       mpTrafo->initArrayCheck(settings->surfaceTrafo);
-       //errMsg("ntlBlenderDumper","mpTrafo inited: "<<(*mpTrafo) );
+       // init world
 }
 
 
@@ -50,7 +41,6 @@ ntlBlenderDumper::ntlBlenderDumper(elbeemSimulationSettings *settings) :
  *****************************************************************************/
 ntlBlenderDumper::~ntlBlenderDumper()
 {
-       delete mpTrafo;
        debMsgStd("ntlBlenderDumper",DM_NOTIFY, "ntlBlenderDumper done", 10);
 }
 
@@ -67,16 +57,6 @@ int ntlBlenderDumper::renderScene( void )
        debugOut = false;
 #endif // ELBEEM_PLUGIN==1
 
-       // output path
-       /*std::ostringstream ecrpath("");
-       ecrpath << "/tmp/ecr_" << getpid() <<"/";
-       // make sure the dir exists
-       std::ostringstream ecrpath_create("");
-       ecrpath_create << "mkdir " << ecrpath.str();
-       system( ecrpath_create.str().c_str() );
-       // */
-
-       vector<string> hideObjs; // geom shaders to hide 
        vector<string> gmName;   // gm names
        vector<string> gmMat;    // materials for gm
        int numGMs = 0;                                  // no. of .obj models created
@@ -91,7 +71,6 @@ int ntlBlenderDumper::renderScene( void )
   vector<ntlTriangle> Triangles;
   vector<ntlVec3Gfx>  Vertices;
   vector<ntlVec3Gfx>  VertNormals;
-       //errMsg("ntlBlenderDumper","mpTrafo : "<<(*mpTrafo) );
 
        /* init geometry array, first all standard objects */
        int idCnt = 0;          // give IDs to objects
@@ -106,14 +85,16 @@ int ntlBlenderDumper::renderScene( void )
                }
                if(tid & GEOCLASSTID_SHADER) {
                        ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter); //dynamic_cast<ntlGeometryShader*>(*iter);
-                       hideObjs.push_back( (*iter)->getName() );
-                       geoshad->notifyShaderOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename());
+                       std::string outname = geoshad->getOutFilename();
+                       if(outname.length()<1) outname = mpGlob->getOutFilename();
+                       geoshad->notifyShaderOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,outname);
+
                        for (vector<ntlGeometryObject*>::iterator siter = geoshad->getObjectsBegin();
                                        siter != geoshad->getObjectsEnd();
                                        siter++) {
                                if(debugOut) debMsgStd("ntlBlenderDumper::BuildScene",DM_MSG,"added shader geometry "<<(*siter)->getName(), 8);
 
-                               (*siter)->notifyOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,glob->getOutFilename(), this->mSimulationTime);
+                               (*siter)->notifyOfDump(DUMP_FULLGEOMETRY, glob->getAniCount(),nrStr,outname, this->mSimulationTime);
                                bool doDump = false;
                                bool isPreview = false;
                                // only dump final&preview surface meshes
@@ -130,8 +111,7 @@ int ntlBlenderDumper::renderScene( void )
                                Vertices.clear();
                                VertNormals.clear();
                                (*siter)->initialize( mpGlob );
-                               //int vstart = mVertNormals.size()-1;
-                               (*siter)->getTriangles(&Triangles, &Vertices, &VertNormals, idCnt);
+                               (*siter)->getTriangles(this->mSimulationTime, &Triangles, &Vertices, &VertNormals, idCnt);
 
                                idCnt ++;
 
@@ -139,8 +119,8 @@ int ntlBlenderDumper::renderScene( void )
 
                                // dump to binary file
                                std::ostringstream boutfilename("");
-                               //boutfilename << ecrpath.str() << glob->getOutFilename() <<"_"<< (*siter)->getName() <<"_" << nrStr << ".obj";
-                               boutfilename << glob->getOutFilename() <<"_"<< (*siter)->getName() <<"_" << nrStr;
+                               //boutfilename << ecrpath.str() << outname <<"_"<< (*siter)->getName() <<"_" << nrStr << ".obj";
+                               boutfilename << outname <<"_"<< (*siter)->getName() <<"_" << nrStr;
                                if(debugOut) debMsgStd("ntlBlenderDumper::renderScene",DM_MSG,"B-Dumping: "<< (*siter)->getName() 
                                                <<", triangles:"<<Triangles.size()<<", vertices:"<<Vertices.size()<<
                                                " to "<<boutfilename.str() , 7);
@@ -164,7 +144,6 @@ int ntlBlenderDumper::renderScene( void )
                                                        // returns smoothed velocity, scaled by frame time
                                                        ntlVec3Gfx v = lbm->getVelocityAt( Vertices[i][0], Vertices[i][1], Vertices[i][2] );
                                                        // translation not necessary, test rotation & scaling?
-                                                       //? v = (*mpTrafo) * v;
                                                        for(int j=0; j<3; j++) {
                                                                float vertp = v[j];
                                                                //if(i<20) errMsg("ntlBlenderDumper","DUMP_VEL final "<<i<<" = "<<v);
@@ -182,10 +161,18 @@ int ntlBlenderDumper::renderScene( void )
                                        errMsg("ntlBlenderDumper::renderScene","Unable to open output '"<<boutfilename<<"' ");
                                        return 1; }
                                
-                               // transform into source space
-                               for(size_t i=0; i<Vertices.size(); i++) {
-                                       Vertices[i] = (*mpTrafo) * Vertices[i];
+                               //! 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];
+                                       }
                                }
+
                                // write to file
                                int numVerts;
                                if(sizeof(numVerts)!=4) { errMsg("ntlBlenderDumper::renderScene","Invalid int size"); return 1; }
index 64f9f81..57f155f 100644 (file)
@@ -9,15 +9,13 @@
 #ifndef NTL_BLENDERDUMPER_H
 #include "ntl_world.h"
 
-template<class Scalar> class ntlMatrix4x4;
-
 class ntlBlenderDumper :
        public ntlWorld
 {
 public:
   /*! Constructor */
+  ntlBlenderDumper();
   ntlBlenderDumper(string filename, bool commandlineMode);
-  ntlBlenderDumper(elbeemSimulationSettings *);
   /*! Destructor */
   virtual ~ntlBlenderDumper( void );
 
@@ -26,8 +24,6 @@ public:
 
 protected:
 
-       //! transform matrix
-       ntlMatrix4x4<gfxReal> *mpTrafo;
 };
 
 #define NTL_BLENDERDUMPER_H
index bcc5a32..0b7a4e9 100644 (file)
@@ -28,8 +28,7 @@ class ntlGeometryClass
                //! Default constructor
                inline ntlGeometryClass() :
                        mVisible( 1 ), mName( "[ObjNameUndef]" ),
-                       mObjectId(-1),
-                       mpAttrs( NULL ) 
+                       mObjectId(-1), mpAttrs( NULL ), mGeoInitId(-1) 
                { 
                                mpAttrs = new AttributeList("objAttrs"); 
                };
@@ -76,6 +75,11 @@ class ntlGeometryClass
                /*! GUI - notify object that mouse was clicked at last pos */
                virtual void setMouseClick() { /* do nothing by default */ }
 
+               /*! Returns the geo init id */
+               inline void setGeoInitId(int set) { mGeoInitId=set; }
+               /*! Returns the geo init id */
+               inline int getGeoInitId() const { return mGeoInitId; }
+
        protected:
 
                /*! Object visible on/off */
@@ -90,6 +94,10 @@ class ntlGeometryClass
                /*! configuration attributes */
                AttributeList *mpAttrs;
 
+               /* fluid init data */
+               /*! id of fluid init (is used in solver initialization), additional data stored only for objects */
+               int mGeoInitId;
+
        private:
 
 };
index 3c2e05a..598abbd 100644 (file)
@@ -26,7 +26,11 @@ ntlGeometryObjModel::ntlGeometryObjModel( void ) :
        ntlGeometryObject(),
        mvStart( 0.0 ), mvEnd( 1.0 ),
        mLoaded( false ),
-       mTriangles(), mVertices(), mNormals()
+       mTriangles(), mVertices(), mNormals(),
+       mcAniVerts(), mcAniNorms(), 
+       mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.),
+       mvCPSStart(-10000.), mvCPSEnd(10000.),
+       mCPSFilename(""), mCPSWidth(0.1), mCPSTimestep(1.)
 {
 }
 
@@ -41,6 +45,13 @@ ntlGeometryObjModel::~ntlGeometryObjModel()
 }
 
 
+/*! is the mesh animated? */
+bool ntlGeometryObjModel::getMeshAnimated() { 
+       const bool ret = (mcAniVerts.getSize()>1); 
+       //errMsg("getMeshAnimated","ret="<<ret<<", size="<<mcAniVerts.getSize() );
+       return ret;
+}
+
 /*****************************************************************************/
 /* Init attributes etc. of this object */
 /*****************************************************************************/
@@ -57,14 +68,14 @@ void ntlGeometryObjModel::initialize(ntlRenderGlobals *glob)
        mFilename = mpAttrs->readString("filename", mFilename,"ntlGeometryObjModel", "mFilename", true);
 
        if(mFilename == "") {
-               errMsg("ntlGeometryObjModel::getTriangles","Filename not given!");
+               errMsg("ntlGeometryObjModel::initialize","Filename not given!");
                return;
        }
 
        const char *suffix = strrchr(mFilename.c_str(), '.');
        if (suffix) {
                if (!strncasecmp(suffix, ".obj", 4)) {
-                       errMsg("ntlGeometryObjModel::getTriangles",".obj files not supported!");
+                       errMsg("ntlGeometryObjModel::initialize",".obj files not supported!");
                        return;
                } else if (!strncasecmp(suffix, ".gz", 3)) { 
                        //mType = 1; // assume its .bobj.gz
@@ -73,18 +84,33 @@ void ntlGeometryObjModel::initialize(ntlRenderGlobals *glob)
                }
        }
 
+       if(getAttributeList()->exists("ani_times") || (!mcAniTimes.isInited()) ) {
+               mcAniTimes = mpAttrs->readChannelFloat("ani_times");
+       }
+       mAniTimeScale = mpAttrs->readFloat("ani_timescale", mAniTimeScale,"ntlGeometryObjModel", "mAniTimeScale", false);
+       mAniTimeOffset = mpAttrs->readFloat("ani_timeoffset", mAniTimeOffset,"ntlGeometryObjModel", "mAniTimeOffset", false);
+
+       mCPSWidth = mpAttrs->readFloat("cps_width", mCPSWidth,"ntlGeometryObjModel", "mCPSWidth", false);
+       mCPSTimestep = mpAttrs->readFloat("cps_timestep", mCPSTimestep,"ntlGeometryObjModel", "mCPSTimestep", false);
+       mvCPSStart = vec2G( mpAttrs->readVec3d("cps_start", vec2D(mvCPSStart),"ntlGeometryObjModel", "mvCPSStart", false));
+       mvCPSEnd = vec2G( mpAttrs->readVec3d("cps_end", vec2D(mvCPSEnd),"ntlGeometryObjModel", "mvCPSEnd", false));
+
        // continue with standard obj
        if(loadBobjModel(mFilename)==0) mLoaded=1;
        if(!mLoaded) {
                debMsgStd("ntlGeometryObjModel",DM_WARNING,"Unable to load object file '"<<mFilename<<"' !", 0);
        }
+       if(getMeshAnimated()) {
+               this->mIsAnimated = true;
+       }
 }
 
 /******************************************************************************
  * init model from given vertex and triangle arrays 
  *****************************************************************************/
 
-int ntlGeometryObjModel::initModel(int numVertices, float *vertices, int numTriangles, int *triangles)
+int ntlGeometryObjModel::initModel(int numVertices, float *vertices, int numTriangles, int *triangles,
+               int channelSize, float *channelVertices)
 {
        mVertices.clear();
        mVertices.resize( numVertices );
@@ -96,10 +122,57 @@ int ntlGeometryObjModel::initModel(int numVertices, float *vertices, int numTria
 
        mTriangles.clear();
        mTriangles.resize( 3*numTriangles );
+       int triangleErrs=0;
        for(int i=0; i<numTriangles; i++) {
-               mTriangles[3*i+0] = triangles[i*3+0];
-               mTriangles[3*i+1] = triangles[i*3+1];
-               mTriangles[3*i+2] = triangles[i*3+2];
+               for(int j=0;j<3;j++) {
+                       mTriangles[3*i+j] = triangles[i*3+j];
+                       if(mTriangles[3*i+j]<0) { mTriangles[3*i+j]=0; triangleErrs++; }
+                       if(mTriangles[3*i+j]>=numVertices) { mTriangles[3*i+j]=0; triangleErrs++; }
+               }
+       }
+       if(triangleErrs>0) {
+               errMsg("ntlGeometryObjModel::initModel","Triangle errors occurred ("<<triangleErrs<<")!");
+       }
+
+       //fprintf(stderr,"initModel DEBUG %d \n",channelSize);
+       debMsgStd("ntlGeometryObjModel::initModel",DM_MSG, "Csize:"<<channelSize<<", Cvert:"<<(long)(channelVertices) ,10);
+       if(channelVertices && (channelSize>0)) {
+               vector<ntlSetVec3f> aniverts;
+               vector<ntlSetVec3f> aninorms;
+               vector<double> anitimes;
+               aniverts.clear();
+               aninorms.clear();
+               anitimes.clear();
+               for(int frame=0; frame<channelSize; frame++) {
+                       ntlSetVec3f averts; averts.mVerts.clear();
+                       ntlSetVec3f anorms; averts.mVerts.clear();
+                       int setsize = (3*numVertices+1);
+
+                       ntlVec3Gfx p(0.),n(1.);
+                       for(int i=0; i<numVertices; i++) {
+                               for(int j=0; j<3; j++) p[j] = channelVertices[frame*setsize+ 3*i +j];
+                               averts.mVerts.push_back(p);
+                               anorms.mVerts.push_back(p);
+                               //debMsgStd("ntlGeometryObjModel::initModel",DM_MSG, "Frame:"<<frame<<",i:"<<i<<" "<<p,10);
+                       }
+                       if( ((int)averts.mVerts.size()==numVertices) && 
+                               ((int)anorms.mVerts.size()==numVertices) ) {
+                               aniverts.push_back(averts);
+                               aninorms.push_back(anorms);
+                               double time = (double)channelVertices[frame*setsize+ setsize-1];
+                               anitimes.push_back(time);
+                       } else {
+                               errMsg("ntlGeometryObjModel::initModel","Invalid mesh, obj="<<this->getName()<<" frame="<<frame<<" verts="<<averts.mVerts.size()<<"/"<<numVertices<<". Skipping...");
+                       }
+                       //debMsgStd("ntlGeometryObjModel::initModel",DM_MSG, "Frame:"<<frame<<" at t="<<time,10);
+               }
+
+               mcAniVerts = AnimChannel<ntlSetVec3f>(aniverts,anitimes);
+               mcAniNorms = AnimChannel<ntlSetVec3f>(aninorms,anitimes);
+               debMsgStd("ntlGeometryObjModel::initModel",DM_MSG, "Ani sets inited: "<< mcAniVerts.accessValues().size() <<","<<mcAniNorms.accessValues().size() <<" ", 1 );
+       }
+       if(getMeshAnimated()) {
+               this->mIsAnimated = true;
        }
 
        // inited, no need to parse attribs etc.
@@ -107,6 +180,48 @@ int ntlGeometryObjModel::initModel(int numVertices, float *vertices, int numTria
        return 0;
 }
 
+/*! init triangle divisions */
+void ntlGeometryObjModel::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri) {
+       // warning - copied from geomobj calc!
+       errMsg("ntlGeometryObjModel","calcTriangleDivs special!");
+       mTriangleDivs1.resize( tris.size() );
+       mTriangleDivs2.resize( tris.size() );
+       mTriangleDivs3.resize( tris.size() );
+       for(size_t i=0; i<tris.size(); i++) {
+               ntlVec3Gfx p0 = verts[ tris[i].getPoints()[0] ];
+               ntlVec3Gfx p1 = verts[ tris[i].getPoints()[1] ];
+               ntlVec3Gfx p2 = verts[ tris[i].getPoints()[2] ];
+               ntlVec3Gfx side1 = p1 - p0;
+               ntlVec3Gfx side2 = p2 - p0;
+               ntlVec3Gfx side3 = p1 - p2;
+               int divs1=0, divs2=0, divs3=0;
+               if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
+               if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
+               //if(normNoSqrt(side3) > fsTri*fsTri) { divs3 = (int)(norm(side3)/fsTri); }
+
+               // special handling
+               // warning, requires objmodel triangle treatment (no verts dups)
+               if(getMeshAnimated()) {
+                       vector<ntlSetVec3f> &sverts = mcAniVerts.accessValues();
+                       for(int s=0; s<(int)sverts.size(); s++) {
+                               p0 = sverts[s].mVerts[ tris[i].getPoints()[0] ];
+                               p1 = sverts[s].mVerts[ tris[i].getPoints()[1] ];
+                               p2 = sverts[s].mVerts[ tris[i].getPoints()[2] ];
+                               side1 = p1 - p0; side2 = p2 - p0; side3 = p1 - p2;
+                               int tdivs1=0, tdivs2=0, tdivs3=0;
+                               if(normNoSqrt(side1) > fsTri*fsTri) { tdivs1 = (int)(norm(side1)/fsTri); }
+                               if(normNoSqrt(side2) > fsTri*fsTri) { tdivs2 = (int)(norm(side2)/fsTri); }
+                               if(tdivs1>divs1) divs1=tdivs1;
+                               if(tdivs2>divs2) divs2=tdivs2;
+                               if(tdivs3>divs3) divs3=tdivs3;
+                       }
+               } // */
+               mTriangleDivs1[i] = divs1;
+               mTriangleDivs2[i] = divs2;
+               mTriangleDivs3[i] = divs3;
+       }
+}
+
 
 /******************************************************************************
  * load model from .obj file
@@ -114,7 +229,13 @@ int ntlGeometryObjModel::initModel(int numVertices, float *vertices, int numTria
 
 int ntlGeometryObjModel::loadBobjModel(string filename)
 {
-       const bool debugPrint=false;
+       bool haveAniSets=false;
+       vector<ntlSetVec3f> aniverts;
+       vector<ntlSetVec3f> aninorms;
+       vector<double> anitimes;
+
+       const bool debugPrint=true;
+       const bool debugPrintFull=true;
        gzFile gzf;
        gzf = gzopen(filename.c_str(), "rb");
        if (!gzf) {
@@ -140,6 +261,7 @@ int ntlGeometryObjModel::loadBobjModel(string filename)
                        gzread(gzf, &(x[j]), sizeof( (x[j]) ) ); 
                }
                mVertices[i] = ntlVec3Gfx(x[0],x[1],x[2]);
+               if(debugPrintFull) errMsg("FULLV"," "<<i<<" "<< mVertices[i] );
        }
        if(debugPrint) errMsg("NV"," "<<numVerts<<" "<< mVertices.size() );
 
@@ -157,6 +279,7 @@ int ntlGeometryObjModel::loadBobjModel(string filename)
                        gzread(gzf, &(n[j]), sizeof( (n[j]) ) ); 
                }
                mNormals[i] = ntlVec3Gfx(n[0],n[1],n[2]);
+               if(debugPrintFull) errMsg("FULLN"," "<<i<<" "<< mNormals[i] );
        }
        if(debugPrint) errMsg("NN"," "<<numVerts<<" "<< mNormals.size() );
 
@@ -180,8 +303,99 @@ int ntlGeometryObjModel::loadBobjModel(string filename)
 
        debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' loaded, #Vertices: "<<mVertices.size()<<", #Normals: "<<mNormals.size()<<", #Triangles: "<<(mTriangles.size()/3)<<" ", 1 );
 
+       // try to load animated mesh
+       aniverts.clear();
+       aninorms.clear();
+       anitimes.clear();
+       while(1) {
+               //ntlVec3Gfx check;
+               float x[3];
+               float frameTime=0.;
+               int bytesRead = 0;
+               int numNorms2=-1, numVerts2=-1;
+               //for(int j=0; j<3; j++) {
+                       //x[j] = 0.;
+                       //bytesRead += gzread(gzf, &(x[j]), sizeof(float) ); 
+               //}
+               //check = ntlVec3Gfx(x[0],x[1],x[2]);
+               //if(debugPrint) errMsg("ANI_NV1"," "<<check<<" "<<" bytes:"<<bytesRead );
+               bytesRead += gzread(gzf, &frameTime, sizeof(frameTime) );
+               //if(bytesRead!=3*sizeof(float)) {
+               if(bytesRead!=sizeof(float)) {
+                       debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' no ani sets. ", 10 );
+                       if(anitimes.size()>0) {
+                               // finally init channels and stop reading file
+                               mcAniVerts = AnimChannel<ntlSetVec3f>(aniverts,anitimes);
+                               mcAniNorms = AnimChannel<ntlSetVec3f>(aninorms,anitimes);
+                       }
+                       goto gzreaddone;
+               }
+               bytesRead += gzread(gzf, &numVerts2, sizeof(numVerts2) );
+               haveAniSets=true;
+               // continue to read new set
+               vector<ntlVec3Gfx> vertset;
+               vector<ntlVec3Gfx> normset;
+               vertset.resize(numVerts);
+               normset.resize(numVerts);
+               //vertset[0] = check;
+               if(debugPrintFull) errMsg("FUL1V"," "<<0<<" "<< vertset[0] );
+
+               for(int i=0; i<numVerts; i++) { // start at one!
+                       for(int j=0; j<3; j++) {
+                               bytesRead += gzread(gzf, &(x[j]), sizeof( (x[j]) ) ); 
+                       }
+                       vertset[i] = ntlVec3Gfx(x[0],x[1],x[2]);
+                       if(debugPrintFull) errMsg("FUL2V"," "<<i<<" "<< vertset[i] );
+               }
+               if(debugPrint) errMsg("ANI_VV"," "<<numVerts<<" "<< vertset.size()<<" bytes:"<<bytesRead );
+
+               bytesRead += gzread(gzf, &numNorms2, sizeof(numNorms2) );
+               for(int i=0; i<numVerts; i++) {
+                       for(int j=0; j<3; j++) {
+                               bytesRead += gzread(gzf, &(x[j]), sizeof( (x[j]) ) ); 
+                       }
+                       normset[i] = ntlVec3Gfx(x[0],x[1],x[2]);
+                       if(debugPrintFull) errMsg("FUL2N"," "<<i<<" "<< normset[i] );
+               }
+               if(debugPrint) errMsg("ANI_NV"," "<<numVerts<<","<<numVerts2<<","<<numNorms2<<","<< normset.size()<<" bytes:"<<bytesRead );
+
+               // set ok
+               if(bytesRead== (int)( (numVerts*2*3+1) *sizeof(float)+2*sizeof(int) ) ) {
+                       if(aniverts.size()==0) {
+                               // TODO, ignore first mesh?
+                               double anitime = (double)(frameTime-1.); // start offset!? anitimes.size();
+                               // get for current frame entry
+                               if(mcAniTimes.getSize()>1)  anitime = mcAniTimes.get(anitime); 
+                               anitime = anitime*mAniTimeScale+mAniTimeOffset;
+
+                               anitimes.push_back( anitime );
+                               aniverts.push_back( ntlSetVec3f(mVertices) );
+                               aninorms.push_back( ntlSetVec3f(mNormals) );
+                               if(debugPrint) errMsg("ANI_NV","new set "<<mVertices.size()<<","<< mNormals.size()<<" time:"<<anitime );
+                       }
+                       double anitime = (double)(frameTime); //anitimes.size();
+                       // get for current frame entry
+                       if(mcAniTimes.getSize()>1)  anitime = mcAniTimes.get(anitime); 
+                       anitime = anitime*mAniTimeScale+mAniTimeOffset;
+
+                       anitimes.push_back( anitime );
+                       aniverts.push_back( ntlSetVec3f(vertset) );
+                       aninorms.push_back( ntlSetVec3f(normset) );
+                       if(debugPrint) errMsg("ANI_NV","new set "<<vertset.size()<<","<< normset.size()<<" time:"<<anitime );
+               } else {
+                       errMsg("ntlGeometryObjModel::loadBobjModel","Malformed ani set! Aborting... ("<<bytesRead<<") ");
+                       goto gzreaddone;
+               }
+       } // anim sets */
+
+gzreaddone:
+
+       if(haveAniSets) { 
+               debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' ani sets loaded: "<< mcAniVerts.accessValues().size() <<","<<mcAniNorms.accessValues().size() <<" ", 1 );
+       }
        gzclose( gzf );
        return 0;
+
 gzreaderror:
        mTriangles.clear();
        mVertices.clear();
@@ -196,24 +410,36 @@ gzreaderror:
  * 
  *****************************************************************************/
 void 
-ntlGeometryObjModel::getTriangles( vector<ntlTriangle> *triangles, 
+ntlGeometryObjModel::getTriangles(double t, vector<ntlTriangle> *triangles, 
                                                                                                                        vector<ntlVec3Gfx> *vertices, 
                                                                                                                        vector<ntlVec3Gfx> *normals, int objectId )
 {
        if(!mLoaded) { // invalid type...
                return;
        }
+       if(mcAniVerts.getSize()>1) { mVertices = mcAniVerts.get(t).mVerts; }
+       if(mcAniNorms.getSize()>1) { mNormals  = mcAniNorms.get(t).mVerts; }
+
+       int startvert = vertices->size();
+       vertices->resize( vertices->size() + mVertices.size() );
+       normals->resize( normals->size() + mVertices.size() );
+       for(int i=0; i<(int)mVertices.size(); i++) {
+               (*vertices)[startvert+i] = mVertices[i];
+               (*normals)[startvert+i] = mNormals[i];
+       }
 
+       triangles->reserve(triangles->size() + mTriangles.size() );
        for(int i=0; i<(int)mTriangles.size(); i+=3) {
                int trip[3];
-               trip[0] = mTriangles[i+0];
-               trip[1] = mTriangles[i+1];
-               trip[2] = mTriangles[i+2];
-
-               sceneAddTriangle( 
-                               mVertices[trip[0]], mVertices[trip[1]], mVertices[trip[2]], 
-                               mNormals[trip[0]], mNormals[trip[1]], mNormals[trip[2]], 
-                               ntlVec3Gfx(0.0), 1 , triangles,vertices,normals ); /* normal unused */
+               trip[0] = startvert+mTriangles[i+0];
+               trip[1] = startvert+mTriangles[i+1];
+               trip[2] = startvert+mTriangles[i+2];
+
+               //sceneAddTriangle( 
+                               //mVertices[trip[0]], mVertices[trip[1]], mVertices[trip[2]], 
+                               //mNormals[trip[0]], mNormals[trip[1]], mNormals[trip[2]], 
+                               //ntlVec3Gfx(0.0), 1 , triangles,vertices,normals ); /* normal unused */
+               sceneAddTriangleNoVert( trip, ntlVec3Gfx(0.0), 1 , triangles ); /* normal unused */
        }
        objectId = -1; // remove warning
        // bobj
index 929b8c5..f911279 100644 (file)
@@ -28,16 +28,24 @@ class ntlGeometryObjModel : public ntlGeometryObject
                /*! Filename setting etc. */
                virtual void initialize(ntlRenderGlobals *glob);
 
+               /*! is the mesh animated? */
+               virtual bool getMeshAnimated();
 
                /* create triangles from obj */
-               virtual void getTriangles( vector<ntlTriangle> *triangles, 
+               virtual void getTriangles(double t,  vector<ntlTriangle> *triangles, 
                                vector<ntlVec3Gfx> *vertices, 
                                vector<ntlVec3Gfx> *normals, int objectId );
 
                /*! load model from .bobj file, returns !=0 upon error */
                int loadBobjModel(string filename);
                /*! init model from given vertex and triangle arrays */
-               int initModel(int numVertices, float *vertices, int numTriangles, int *triangles);
+               int initModel(int numVertices, float *vertices, int numTriangles, int *triangles,
+                               int channelSize, float *channelVertices);
+               /*! init triangle divisions */
+               virtual void calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri);
+
+               /*! do ani mesh CPS */
+               void calculateCPS(string filename);
 
        private:
 
@@ -55,6 +63,20 @@ class ntlGeometryObjModel : public ntlGeometryObject
                vector<ntlVec3Gfx> mVertices;
                vector<ntlVec3Gfx> mNormals;
 
+               /*! animated channels for vertices, if given will override getris by default */
+               AnimChannel<ntlSetVec3f> mcAniVerts;
+               AnimChannel<ntlSetVec3f> mcAniNorms;
+               /*! map entrie of anim mesh to sim times */
+               AnimChannel<double> mcAniTimes;
+               /*! timing mapping & offset for config files */
+               double mAniTimeScale, mAniTimeOffset;
+
+               /*! ani mesh cps params */
+               ntlVec3Gfx mvCPSStart, mvCPSEnd;
+               string mCPSFilename;
+               gfxReal mCPSWidth, mCPSTimestep;
+
+
        public:
 
                /* Access methods */
index 5ab13a9..90c3b34 100644 (file)
 ntlGeometryObject::ntlGeometryObject() :
        mIsInitialized(false), mpMaterial( NULL ),
        mMaterialName( "default" ),
-       mCastShadows( 1 ),
-       mReceiveShadows( 1 ),
-       mGeoInitId( -1 ), mGeoInitType( 0 ), 
+       mCastShadows( 1 ), mReceiveShadows( 1 ),
+       mGeoInitType( 0 ), 
        mInitialVelocity(0.0), mcInitialVelocity(0.0), mLocalCoordInivel(false),
        mGeoInitIntersect(false),
        mGeoPartSlipValue(0.0),
-       mOnlyThinInit(false),
+       mVolumeInit(VOLUMEINIT_VOLUME),
        mInitialPos(0.),
        mcTrans(0.), mcRot(0.), mcScale(1.),
        mIsAnimated(false),
-       mMovPoints(), mMovNormals(),
+       mMovPoints(), //mMovNormals(),
        mHaveCachedMov(false),
-       mCachedMovPoints(), mCachedMovNormals(),
+       mCachedMovPoints(), //mCachedMovNormals(),
+       mTriangleDivs1(), mTriangleDivs2(), mTriangleDivs3(),
        mMovPntsInited(-100.0), mMaxMovPnt(-1),
        mcGeoActive(1.)
 { 
@@ -50,6 +50,30 @@ ntlGeometryObject::~ntlGeometryObject()
 {
 }
 
+/*! is the mesh animated? */
+bool ntlGeometryObject::getMeshAnimated() {
+       // off by default, on for e.g. ntlGeometryObjModel
+       return false; 
+}
+
+/*! init object anim flag */
+bool ntlGeometryObject::checkIsAnimated() {
+       if(    (mcTrans.accessValues().size()>1)  // VALIDATE
+           || (mcRot.accessValues().size()>1) 
+           || (mcScale.accessValues().size()>1) 
+           || (mcGeoActive.accessValues().size()>1) 
+           || (mcInitialVelocity.accessValues().size()>1) 
+               ) {
+               mIsAnimated = true;
+       }
+
+       // fluid objects always have static init!
+       if(mGeoInitType==FGI_FLUID) {
+               mIsAnimated=false;
+       }
+       return mIsAnimated;
+}
+
 /*****************************************************************************/
 /* Init attributes etc. of this object */
 /*****************************************************************************/
@@ -77,10 +101,10 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
        // init material, always necessary
        searchMaterial( glob->getMaterials() );
        
-       mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"ntlGeometryObject", "mGeoInitId", false);
+       this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"ntlGeometryObject", "mGeoInitId", false);
        mGeoInitIntersect = mpAttrs->readInt("geoinit_intersect", mGeoInitIntersect,"ntlGeometryObject", "mGeoInitIntersect", false);
        string ginitStr = mpAttrs->readString("geoinittype", "", "ntlGeometryObject", "mGeoInitType", false);
-       if(mGeoInitId>=0) {
+       if(this->mGeoInitId>=0) {
                bool gotit = false;
                for(int i=0; i<GEOINIT_STRINGS; i++) {
                        if(ginitStr== initStringStrs[i]) {
@@ -95,10 +119,10 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
                }
        }
 
-       int geoActive = mpAttrs->readInt("geoinitactive", 1,"ntlGeometryObject", "mGeoInitId", false);
+       int geoActive = mpAttrs->readInt("geoinitactive", 1,"ntlGeometryObject", "geoActive", false);
        if(!geoActive) {
                // disable geo init again...
-               mGeoInitId = -1;
+               this->mGeoInitId = -1;
        }
        mInitialVelocity  = vec2G( mpAttrs->readVec3d("initial_velocity", vec2D(mInitialVelocity),"ntlGeometryObject", "mInitialVelocity", false));
        if(getAttributeList()->exists("initial_velocity") || (!mcInitialVelocity.isInited()) ) {
@@ -109,7 +133,11 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
        mLocalCoordInivel = mpAttrs->readBool("geoinit_localinivel", mLocalCoordInivel,"ntlGeometryObject", "mLocalCoordInivel", false);
 
        mGeoPartSlipValue = mpAttrs->readFloat("geoinit_partslip", mGeoPartSlipValue,"ntlGeometryObject", "mGeoPartSlipValue", false);
+       bool mOnlyThinInit = false; // deprecated!
        mOnlyThinInit     = mpAttrs->readBool("geoinit_onlythin", mOnlyThinInit,"ntlGeometryObject", "mOnlyThinInit", false);
+       if(mOnlyThinInit) mVolumeInit = VOLUMEINIT_SHELL;
+       mVolumeInit     = mpAttrs->readInt("geoinit_volumeinit", mVolumeInit,"ntlGeometryObject", "mVolumeInit", false);
+       if((mVolumeInit<VOLUMEINIT_VOLUME)||(mVolumeInit>VOLUMEINIT_BOTH)) mVolumeInit = VOLUMEINIT_VOLUME;
 
        // override cfg types
        mVisible = mpAttrs->readBool("visible", mVisible,"ntlGeometryObject", "mVisible", false);
@@ -141,17 +169,18 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
        // always use channel
        if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel<double>(geoactive); }
 
-       if(    (mcTrans.accessValues().size()>1)  // VALIDATE
-           || (mcRot.accessValues().size()>1) 
-           || (mcScale.accessValues().size()>1) 
-           || (mcGeoActive.accessValues().size()>1) 
-           || (mcInitialVelocity.accessValues().size()>1) 
-               ) {
-               mIsAnimated = true;
-       }
+       //if(    (mcTrans.accessValues().size()>1)  // VALIDATE
+           //|| (mcRot.accessValues().size()>1) 
+           //|| (mcScale.accessValues().size()>1) 
+           //|| (mcGeoActive.accessValues().size()>1) 
+           //|| (mcInitialVelocity.accessValues().size()>1) 
+               //) {
+               //mIsAnimated = true;
+       //}
+       checkIsAnimated();
 
        mIsInitialized = true;
-       debMsgStd("ntlGeometryObject::initialize",DM_MSG,"GeoObj '"<<this->getName()<<"': visible="<<this->mVisible<<" gid="<<mGeoInitId<<" gtype="<<mGeoInitType<<","<<ginitStr<<
+       debMsgStd("ntlGeometryObject::initialize",DM_MSG,"GeoObj '"<<this->getName()<<"': visible="<<this->mVisible<<" gid="<<this->mGeoInitId<<" gtype="<<mGeoInitType<<","<<ginitStr<<
                        " gvel="<<mInitialVelocity<<" gisect="<<mGeoInitIntersect, 10); // debug
 }
 
@@ -222,16 +251,11 @@ void ntlGeometryObject::sceneAddTriangle(
                if(getVisible()){ flag |= TRI_GEOMETRY; }
                if(getCastShadows() ) { 
                        flag |= TRI_CASTSHADOWS; } 
-               if( (getMaterial()->getMirror()>0.0) ||  
-                               (getMaterial()->getTransparence()>0.0) ||  
-                               (getMaterial()->getFresnel()>0.0) ) { 
-                       flag |= TRI_MAKECAUSTICS; } 
-               else { 
-                       flag |= TRI_NOCAUSTICS; } 
                
                /* init geo init id */
                int geoiId = getGeoInitId(); 
-               if((geoiId > 0) && (!mOnlyThinInit) && (!mIsAnimated)) { 
+               //if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME) && (!mIsAnimated)) { 
+               if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME)) { 
                        flag |= (1<< (geoiId+4)); 
                        flag |= mGeoInitType; 
                } 
@@ -246,6 +270,38 @@ void ntlGeometryObject::sceneAddTriangle(
        } /* normals check*/ 
 }
 
+void ntlGeometryObject::sceneAddTriangleNoVert(int *trips,
+               ntlVec3Gfx trin, bool smooth,
+               vector<ntlTriangle> *triangles) {
+       ntlTriangle tri;
+               
+       tri.getPoints()[0] = trips[0];
+       tri.getPoints()[1] = trips[1];
+       tri.getPoints()[2] = trips[2];
+
+       // same as normal sceneAddTriangle
+
+       /* init flags from ntl_ray.h */
+       int flag = 0; 
+       if(getVisible()){ flag |= TRI_GEOMETRY; }
+       if(getCastShadows() ) { 
+               flag |= TRI_CASTSHADOWS; } 
+
+       /* init geo init id */
+       int geoiId = getGeoInitId(); 
+       if((geoiId > 0) && (mVolumeInit&VOLUMEINIT_VOLUME)) { 
+               flag |= (1<< (geoiId+4)); 
+               flag |= mGeoInitType; 
+       } 
+       /*errMsg("ntlScene::addTriangle","DEBUG flag="<<convertFlags2String(flag) ); */ 
+       tri.setFlags( flag );
+
+       /* triangle normal missing */
+       tri.setNormal( trin );
+       tri.setSmoothNormals( smooth );
+       tri.setObjectId( this->mObjectId );
+       triangles->push_back( tri ); 
+}
 
 
 /******************************************************************************/
@@ -284,14 +340,15 @@ void ntlGeometryObject::initChannels(
        if((act)&&(nAct>0)) {      ADD_CHANNEL_FLOAT(mcGeoActive, nAct, act); }
        if((ivel)&&(nIvel>0)) {    ADD_CHANNEL_VEC(mcInitialVelocity, nIvel, ivel); }
 
-       if(    (mcTrans.accessValues().size()>1)  // VALIDATE
-           || (mcRot.accessValues().size()>1) 
-           || (mcScale.accessValues().size()>1) 
-           || (mcGeoActive.accessValues().size()>1) 
-           || (mcInitialVelocity.accessValues().size()>1) 
-               ) {
-               mIsAnimated = true;
-       }
+       //if(    (mcTrans.accessValues().size()>1)  // VALIDATE
+           //|| (mcRot.accessValues().size()>1) 
+           //|| (mcScale.accessValues().size()>1) 
+           //|| (mcGeoActive.accessValues().size()>1) 
+           //|| (mcInitialVelocity.accessValues().size()>1) 
+               //) {
+               //mIsAnimated = true;
+       //}
+       checkIsAnimated();
        if(debugInitc) { 
                debMsgStd("ntlGeometryObject::initChannels",DM_MSG,getName()<<
                                " nt:"<<mcTrans.accessValues().size()<<" nr:"<<mcRot.accessValues().size()<<
@@ -336,7 +393,7 @@ void ntlGeometryObject::applyTransformation(double t, vector<ntlVec3Gfx> *verts,
            || (!mHaveCachedMov)
                ) {
                // transformation is animated, continue
-               ntlVec3Gfx pos = mcTrans.get(t);
+               ntlVec3Gfx pos = getTranslation(t); 
                ntlVec3Gfx scale = mcScale.get(t);
                ntlVec3Gfx rot = mcRot.get(t);
                ntlMat4Gfx rotMat;
@@ -347,6 +404,7 @@ void ntlGeometryObject::applyTransformation(double t, vector<ntlVec3Gfx> *verts,
                        (*verts)[i] *= scale;
                        (*verts)[i] = rotMat * (*verts)[i];
                        (*verts)[i] += pos;
+                       //if(i<10) errMsg("ntlGeometryObject::applyTransformation"," v"<<i<<"/"<<vend<<"="<<(*verts)[i]);
                }
                if(norms) {
                        for(int i=vstart; i<vend; i++) {
@@ -355,30 +413,62 @@ 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 ");
+       }
+}
+
+/*! init triangle divisions */
+void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri) {
+       mTriangleDivs1.resize( tris.size() );
+       mTriangleDivs2.resize( tris.size() );
+       mTriangleDivs3.resize( tris.size() );
+       for(size_t i=0; i<tris.size(); i++) {
+               const ntlVec3Gfx p0 = verts[ tris[i].getPoints()[0] ];
+               const ntlVec3Gfx p1 = verts[ tris[i].getPoints()[1] ];
+               const ntlVec3Gfx p2 = verts[ tris[i].getPoints()[2] ];
+               const ntlVec3Gfx side1 = p1 - p0;
+               const ntlVec3Gfx side2 = p2 - p0;
+               const ntlVec3Gfx side3 = p1 - p2;
+               int divs1=0, divs2=0, divs3=0;
+               if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
+               if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
+               //if(normNoSqrt(side3) > fsTri*fsTri) { divs3 = (int)(norm(side3)/fsTri); }
+               /*if(getMeshAnimated()) {
+                       vector<ntlSetVec3f> *verts =mcAniVerts.accessValues();
+                       for(int s=0; s<verts->size(); s++) {
+                               int tdivs1=0, tdivs2=0, tdivs3=0;
+                               if(normNoSqrt(side1) > fsTri*fsTri) { tdivs1 = (int)(norm(side1)/fsTri); }
+                               if(normNoSqrt(side2) > fsTri*fsTri) { tdivs2 = (int)(norm(side2)/fsTri); }
+                               if(tdivs1>divs1) divs1=tdivs1;
+                               if(tdivs2>divs2) divs2=tdivs2;
+                       }
+               }*/
+               mTriangleDivs1[i] = divs1;
+               mTriangleDivs2[i] = divs2;
+               mTriangleDivs3[i] = divs3;
        }
 }
 
 /*! Prepare points for moving objects */
-void ntlGeometryObject::initMovingPoints(gfxReal featureSize) {
-       if(mMovPntsInited==featureSize) return;
+void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
+       if((mMovPntsInited==featureSize)&&(!getMeshAnimated())) return;
        const bool debugMoinit=false;
 
+       //vector<ntlVec3Gfx> movNormals;
        vector<ntlTriangle> triangles; 
        vector<ntlVec3Gfx> vertices; 
-       vector<ntlVec3Gfx> normals; 
+       vector<ntlVec3Gfx> vnormals; 
        int objectId = 1;
-       this->getTriangles(&triangles,&vertices,&normals,objectId);
+       this->getTriangles(time, &triangles,&vertices,&vnormals,objectId);
        
        mMovPoints.clear(); //= vertices;
-       mMovNormals.clear(); //= normals;
+       //movNormals.clear(); //= vnormals;
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" has v:"<<vertices.size()<<" t:"<<triangles.size() );
        // no points?
        if(vertices.size()<1) {
                mMaxMovPnt=-1;
                return; 
        }
-
        ntlVec3f maxscale = channelFindMaxVf(mcScale);
        float maxpart = ABS(maxscale[0]);
        if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
@@ -388,8 +478,12 @@ void ntlGeometryObject::initMovingPoints(gfxReal featureSize) {
        const gfxReal fsTri = featureSize*0.5 *scaleFac;
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","maxscale:"<<maxpart<<" featureSize:"<<featureSize<<" fsTri:"<<fsTri );
 
+       if(mTriangleDivs1.size()!=triangles.size()) {
+               calcTriangleDivs(vertices,triangles,fsTri);
+       }
+
        // debug: count points to init
-       if(debugMoinit) {
+       /*if(debugMoinit) {
                errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" estimating...");
                int countp=vertices.size()*2;
                for(size_t i=0; i<triangles.size(); i++) {
@@ -412,7 +506,7 @@ void ntlGeometryObject::initMovingPoints(gfxReal featureSize) {
                        }
                }
                errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" requires:"<<countp*2);
-       }
+       } // */
 
        bool discardInflowBack = false;
        if( (mGeoInitType==FGI_MBNDINFLOW) && (mcInitialVelocity.accessValues().size()<1) ) discardInflowBack = true;
@@ -422,79 +516,87 @@ void ntlGeometryObject::initMovingPoints(gfxReal featureSize) {
        // init std points
        for(size_t i=0; i<vertices.size(); i++) {
                ntlVec3Gfx p = vertices[ i ];
-               ntlVec3Gfx n = normals[ i ];
+               ntlVec3Gfx n = vnormals[ i ];
                // discard inflow backsides
                //if( (mGeoInitType==FGI_MBNDINFLOW) && (!mIsAnimated)) {
                if(discardInflowBack) { //if( (mGeoInitType==FGI_MBNDINFLOW) && (!mIsAnimated)) {
                        if(dot(mInitialVelocity,n)<0.0) continue;
                }
                mMovPoints.push_back(p);
-               mMovNormals.push_back(n);
+               //movNormals.push_back(n);
+               //errMsg("ntlGeometryObject::initMovingPoints","std"<<i<<" p"<<p<<" n"<<n<<" ");
        }
        // init points & refine...
        for(size_t i=0; i<triangles.size(); i++) {
-               ntlVec3Gfx p0 = vertices[ triangles[i].getPoints()[0] ];
-               ntlVec3Gfx side1 = vertices[ triangles[i].getPoints()[1] ] - p0;
-               ntlVec3Gfx side2 = vertices[ triangles[i].getPoints()[2] ] - p0;
-               int divs1=0, divs2=0;
-               if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
-               if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
-               /* if( (i!=6) &&
-                               (i!=6) ) { divs1=divs2=0; } // DEBUG */
+               int *trips = triangles[i].getPoints();
+               const ntlVec3Gfx p0 = vertices[ trips[0] ];
+               const ntlVec3Gfx side1 = vertices[ trips[1] ] - p0;
+               const ntlVec3Gfx side2 = vertices[ trips[2] ] - p0;
+               int divs1=mTriangleDivs1[i], divs2=mTriangleDivs2[i];
+               //if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
+               //if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
+               const ntlVec3Gfx trinorm = getNormalized(cross(side1,side2))*0.5*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(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);
                                        if(uf+vf>1.0) continue;
-                                       ntlVec3Gfx p = vertices[ triangles[i].getPoints()[0] ] * (1.0-uf-vf)+
-                                               vertices[ triangles[i].getPoints()[1] ]*uf +
-                                               vertices[ triangles[i].getPoints()[2] ]*vf;
-                                       ntlVec3Gfx n = normals[ triangles[i].getPoints()[0] ] * (1.0-uf-vf)+
-                                               normals[ triangles[i].getPoints()[1] ]*uf +
-                                               normals[ triangles[i].getPoints()[2] ]*vf;
-                                       normalize(n);
-                                       //if(mGeoInitType==FGI_MBNDINFLOW) {
+                                       ntlVec3Gfx p = 
+                                               vertices[ trips[0] ] * (1.0-uf-vf)+
+                                               vertices[ trips[1] ] * uf +
+                                               vertices[ trips[2] ] * vf;
+                                       //ntlVec3Gfx n = vnormals[ 
+                                               //trips[0] ] * (1.0-uf-vf)+
+                                               //vnormals[ trips[1] ]*uf +
+                                               //vnormals[ trips[2] ]*vf;
+                                       //normalize(n);
                                        // discard inflow backsides
-                                       if(discardInflowBack) { //if( (mGeoInitType==FGI_MBNDINFLOW) && (!mIsAnimated)) {
-                                               if(dot(mInitialVelocity,n)<0.0) continue;
-                                       }
 
                                        mMovPoints.push_back(p);
-                                       mMovNormals.push_back(n);
+                                       mMovPoints.push_back(p - trinorm);
+                                       //movNormals.push_back(n);
+                                       //movNormals.push_back(trinorm);
+                                       //errMsg("TRINORM","p"<<p<<" n"<<n<<" trin"<<trinorm);
                                }
                        }
                }
        }
 
        // duplicate insides
-       size_t mpsize = mMovPoints.size();
-       for(size_t i=0; i<mpsize; i++) {
-               //normalize(normals[i]);
+       //size_t mpsize = mMovPoints.size();
+       //for(size_t i=0; i<mpsize; i++) {
+               //normalize(vnormals[i]);
                //errMsg("TTAT"," moved:"<<(mMovPoints[i] - mMovPoints[i]*featureSize)<<" org"<<mMovPoints[i]<<" norm"<<mMovPoints[i]<<" fs"<<featureSize);
-               mMovPoints.push_back(mMovPoints[i] - mMovNormals[i]*0.5*featureSize);
-               mMovNormals.push_back(mMovNormals[i]);
-       }
+               //mMovPoints.push_back(mMovPoints[i] - movNormals[i]*0.5*featureSize);
+               //movNormals.push_back(movNormals[i]);
+       //}
 
        // find max point
        mMaxMovPnt = 0;
        gfxReal dist = normNoSqrt(mMovPoints[0]);
-       for(size_t i=0; i<mpsize; i++) {
+       for(size_t i=0; i<mMovPoints.size(); i++) {
                if(normNoSqrt(mMovPoints[i])>dist) {
                        mMaxMovPnt = i;
                        dist = normNoSqrt(mMovPoints[0]);
                }
        }
 
-       if(    (mcTrans.accessValues().size()>1)  // VALIDATE
+       if(    (this-getMeshAnimated())
+      || (mcTrans.accessValues().size()>1)  // VALIDATE
            || (mcRot.accessValues().size()>1) 
            || (mcScale.accessValues().size()>1) 
                ) {
                // also do trafo...
        } else {
                mCachedMovPoints = mMovPoints;
-               mCachedMovNormals = mMovNormals;
-               applyTransformation(0., &mCachedMovPoints, &mCachedMovNormals, 0, mCachedMovPoints.size(), true);
+               //mCachedMovNormals = movNormals;
+               //applyTransformation(time, &mCachedMovPoints, &mCachedMovNormals, 0, mCachedMovPoints.size(), true);
+               applyTransformation(time, &mCachedMovPoints, NULL, 0, mCachedMovPoints.size(), true);
                mHaveCachedMov = true;
                debMsgStd("ntlGeometryObject::initMovingPoints",DM_MSG,"Object "<<getName()<<" cached points ", 7);
        }
@@ -502,18 +604,186 @@ void ntlGeometryObject::initMovingPoints(gfxReal featureSize) {
        mMovPntsInited = featureSize;
        debMsgStd("ntlGeometryObject::initMovingPoints",DM_MSG,"Object "<<getName()<<" inited v:"<<vertices.size()<<"->"<<mMovPoints.size() , 5);
 }
+/*! Prepare points for animated objects,
+ * init both sets, never used cached points 
+ * discardInflowBack ignore */
+void ntlGeometryObject::initMovingPointsAnim(
+                double srctime, vector<ntlVec3Gfx> &srcmovPoints,
+                double dsttime, vector<ntlVec3Gfx> &dstmovPoints,
+                gfxReal featureSize,
+                ntlVec3Gfx geostart, ntlVec3Gfx geoend
+                ) {
+       const bool debugMoinit=false;
+
+       //vector<ntlVec3Gfx> srcmovNormals;
+       //vector<ntlVec3Gfx> dstmovNormals;
+
+       vector<ntlTriangle> srctriangles; 
+       vector<ntlVec3Gfx> srcvertices; 
+       vector<ntlVec3Gfx> unused_normals; 
+       vector<ntlTriangle> dsttriangles; 
+       vector<ntlVec3Gfx> dstvertices; 
+       //vector<ntlVec3Gfx> dstnormals; 
+       int objectId = 1;
+       // TODO optimize? , get rid of normals?
+       unused_normals.clear();
+       this->getTriangles(srctime, &srctriangles,&srcvertices,&unused_normals,objectId);
+       unused_normals.clear();
+       this->getTriangles(dsttime, &dsttriangles,&dstvertices,&unused_normals,objectId);
+       
+       srcmovPoints.clear();
+       dstmovPoints.clear();
+       //srcmovNormals.clear();
+       //dstmovNormals.clear();
+       if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has srcv:"<<srcvertices.size()<<" srct:"<<srctriangles.size() );
+       if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has dstv:"<<dstvertices.size()<<" dstt:"<<dsttriangles.size() );
+       // no points?
+       if(srcvertices.size()<1) {
+               mMaxMovPnt=-1;
+               return; 
+       }
+       if((srctriangles.size() != dsttriangles.size()) ||
+          (srcvertices.size() != dstvertices.size()) ) {
+               errMsg("ntlGeometryObject::initMovingPointsAnim","Invalid triangle numbers! Aborting...");
+               return;
+       }
+       ntlVec3f maxscale = channelFindMaxVf(mcScale);
+       float maxpart = ABS(maxscale[0]);
+       if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
+       if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
+       float scaleFac = 1.0/(maxpart);
+       // TODO - better reinit from time to time?
+       const gfxReal fsTri = featureSize*0.5 *scaleFac;
+       if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","maxscale:"<<maxpart<<" featureSize:"<<featureSize<<" fsTri:"<<fsTri );
+
+       if(mTriangleDivs1.size()!=srctriangles.size()) {
+               calcTriangleDivs(srcvertices,srctriangles,fsTri);
+       }
+
+
+       // init std points
+       for(size_t i=0; i<srcvertices.size(); i++) {
+               srcmovPoints.push_back(srcvertices[i]);
+               //srcmovNormals.push_back(srcnormals[i]);
+       }
+       for(size_t i=0; i<dstvertices.size(); i++) {
+               dstmovPoints.push_back(dstvertices[i]);
+               //dstmovNormals.push_back(dstnormals[i]);
+       }
+       if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<" dst:"<<dstmovPoints.size()<<" " );
+       // init points & refine...
+       for(size_t i=0; i<srctriangles.size(); i++) {
+               const int divs1=mTriangleDivs1[i];
+               const int       divs2=mTriangleDivs2[i];
+               if(divs1+divs2 > 0) {
+                       int *srctrips = srctriangles[i].getPoints();
+                       int *dsttrips = dsttriangles[i].getPoints();
+                       const ntlVec3Gfx srcp0 =    srcvertices[ srctrips[0] ];
+                       const ntlVec3Gfx srcside1 = srcvertices[ srctrips[1] ] - srcp0;
+                       const ntlVec3Gfx srcside2 = srcvertices[ srctrips[2] ] - srcp0;
+                       const ntlVec3Gfx dstp0 =    dstvertices[ dsttrips[0] ];
+                       const ntlVec3Gfx dstside1 = dstvertices[ dsttrips[1] ] - dstp0;
+                       const ntlVec3Gfx dstside2 = dstvertices[ dsttrips[2] ] - dstp0;
+                       const ntlVec3Gfx src_trinorm = getNormalized(cross(srcside1,srcside2))*0.5*featureSize;
+                       const ntlVec3Gfx dst_trinorm = getNormalized(cross(dstside1,dstside2))*0.5*featureSize;
+                       //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);
+                                       if(uf+vf>1.0) continue;
+                                       ntlVec3Gfx srcp = 
+                                               srcvertices[ srctrips[0] ] * (1.0-uf-vf)+
+                                               srcvertices[ srctrips[1] ] * uf +
+                                               srcvertices[ srctrips[2] ] * vf;
+                                       ntlVec3Gfx dstp = 
+                                               dstvertices[ dsttrips[0] ] * (1.0-uf-vf)+
+                                               dstvertices[ dsttrips[1] ] * uf +
+                                               dstvertices[ dsttrips[2] ] * vf;
+
+                                       // cutoffDomain
+                                       if((srcp[0]<geostart[0]) && (dstp[0]<geostart[0])) continue;
+                                       if((srcp[1]<geostart[1]) && (dstp[1]<geostart[1])) continue;
+                                       if((srcp[2]<geostart[2]) && (dstp[2]<geostart[2])) continue;
+                                       if((srcp[0]>geoend[0]  ) && (dstp[0]>geoend[0]  )) continue;
+                                       if((srcp[1]>geoend[1]  ) && (dstp[1]>geoend[1]  )) continue;
+                                       if((srcp[2]>geoend[2]  ) && (dstp[2]>geoend[2]  )) continue;
+                                       
+                                       //ntlVec3Gfx srcn = 
+                                               //srcnormals[ srctriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
+                                               //srcnormals[ srctriangles[i].getPoints()[1] ] * uf +
+                                               //srcnormals[ srctriangles[i].getPoints()[2] ] * vf;
+                                       //normalize(srcn);
+                                       srcmovPoints.push_back(srcp);
+                                       srcmovPoints.push_back(srcp-src_trinorm);
+                                       //srcmovNormals.push_back(srcn);
+
+                                       //ntlVec3Gfx dstn = 
+                                               //dstnormals[ dsttriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
+                                               //dstnormals[ dsttriangles[i].getPoints()[1] ] * uf +
+                                               //dstnormals[ dsttriangles[i].getPoints()[2] ] * vf;
+                                       //normalize(dstn);
+                                       dstmovPoints.push_back(dstp);
+                                       dstmovPoints.push_back(dstp-dst_trinorm);
+                                       //dstmovNormals.push_back(dstn);
+
+                                       /*if(debugMoinit && (i>=0)) errMsg("ntlGeometryObject::initMovingPointsAnim"," "<<
+                                       //srcmovPoints[ srcmovPoints.size()-1 ]<<","<<
+                                       //srcmovNormals[ srcmovNormals.size()-1 ]<<"   "<<
+                                       //dstmovPoints[ dstmovPoints.size()-1 ]<<","<<
+                                       //dstmovNormals[ dstmovNormals.size()-1 ]<<" " 
+                                       (srcmovNormals[ srcmovPoints.size()-1 ]-
+                                       dstmovNormals[ dstmovPoints.size()-1 ])<<" "
+                                       );
+                                       // */
+                               }
+                       }
+               }
+       }
+
+       /*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
+       // duplicate insides
+       size_t mpsize = srcmovPoints.size();
+       for(size_t i=0; i<mpsize; i++) {
+               srcmovPoints.push_back(srcmovPoints[i] - srcmovNormals[i]*1.0*featureSize);
+               //? srcnormals.push_back(srcnormals[i]);
+       }
+       mpsize = dstmovPoints.size();
+       for(size_t i=0; i<mpsize; i++) {
+               dstmovPoints.push_back(dstmovPoints[i] - dstmovNormals[i]*1.0*featureSize);
+               //? dstnormals.push_back(dstnormals[i]);
+       }
+       // */
+
+       /*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
+       for(size_t i=0; i<srcmovPoints.size(); i++) {
+               ntlVec3Gfx p1 = srcmovPoints[i];
+               ntlVec3Gfx p2 = dstmovPoints[i];
+         gfxReal len = norm(p1-p2);
+               if(len>0.01) { errMsg("ntlGeometryObject::initMovingPointsAnim"," i"<<i<<" "<< p1<<" "<<p2<<" "<<len); }
+       } // */
+
+       // find max point not necessary
+       debMsgStd("ntlGeometryObject::initMovingPointsAnim",DM_MSG,"Object "<<getName()<<" inited v:"<<srcvertices.size()<<"->"<<srcmovPoints.size()<<","<<dstmovPoints.size() , 5);
+}
 
 /*! Prepare points for moving objects */
 void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3Gfx> *norms) {
        if(mHaveCachedMov) {
                ret = mCachedMovPoints;
-               if(norms) { *norms = mCachedMovNormals; }
-               errMsg ("ntlGeometryObject::getMovingPoints","Object "<<getName()<<" used cached points ");
+               if(norms) { 
+                       //*norms = mCachedMovNormals; 
+                       errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
+               }
+               //errMsg ("ntlGeometryObject::getMovingPoints","Object "<<getName()<<" used cached points "); // DEBUG
                return;
        }
 
        ret = mMovPoints;
-       if(norms) { *norms = mMovNormals; }
+       if(norms) { 
+               errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
+               //*norms = mMovNormals; 
+       }
 }
 
 
@@ -536,6 +806,31 @@ ntlVec3Gfx ntlGeometryObject::calculateMaxVel(double t1, double t2) {
 /*! get translation at time t*/
 ntlVec3Gfx ntlGeometryObject::getTranslation(double t) {
        ntlVec3Gfx pos = mcTrans.get(t);
+  // DEBUG CP_FORCECIRCLEINIT 1
+       /*if( 
+                       (mName.compare(string("0__ts1"))==0) ||
+                       (mName.compare(string("1__ts1"))==0) ||
+                       (mName.compare(string("2__ts1"))==0) ||
+                       (mName.compare(string("3__ts1"))==0) ||
+                       (mName.compare(string("4__ts1"))==0) ||
+                       (mName.compare(string("5__ts1"))==0) ||
+                       (mName.compare(string("6__ts1"))==0) ||
+                       (mName.compare(string("7__ts1"))==0) ||
+                       (mName.compare(string("8__ts1"))==0) ||
+                       (mName.compare(string("9__ts1"))==0) 
+                       ) { int j=mName[0]-'0';
+                       ntlVec3Gfx ppos(0.); { // DEBUG
+                       const float tscale=10.;
+                       const float tprevo = 0.33;
+                       const ntlVec3Gfx   toff(50,50,0);
+                       const ntlVec3Gfx oscale(30,30,0);
+                       ppos[0] =  cos(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[0] + toff[0];
+                       ppos[1] = -sin(tscale* t - tprevo*(float)j + M_PI -0.1) * oscale[1] + toff[1];
+                       ppos[2] =                               toff[2]; } // DEBUG
+                       pos = ppos;
+                       pos[2] = 0.15;
+       }
+  // DEBUG CP_FORCECIRCLEINIT 1 */
        return pos;
 }
 /*! get active flag time t*/
index ae81d74..5fe9946 100644 (file)
@@ -19,6 +19,10 @@ class ntlTriangle;
 #define DUMP_FULLGEOMETRY 1
 #define DUMP_PARTIAL      2
 
+#define VOLUMEINIT_VOLUME 1
+#define VOLUMEINIT_SHELL  2
+#define VOLUMEINIT_BOTH   (VOLUMEINIT_SHELL|VOLUMEINIT_VOLUME)
+
 class ntlGeometryObject : public ntlGeometryClass
 {
 
@@ -35,7 +39,7 @@ class ntlGeometryObject : public ntlGeometryClass
                virtual void initialize(ntlRenderGlobals *glob);
 
                /*! Get the triangles from this object */
-               virtual void getTriangles( vector<ntlTriangle> *triangles, 
+               virtual void getTriangles(double t, vector<ntlTriangle> *triangles, 
                                vector<ntlVec3Gfx> *vertices, 
                                vector<ntlVec3Gfx> *normals, int objectId ) = 0;
                
@@ -65,12 +69,8 @@ class ntlGeometryObject : public ntlGeometryClass
                /*! Returns the cast shadows attribute */
                inline int getCastShadows() const { return mCastShadows; }
 
-               /*! Returns the geo init id */
-               inline void setGeoInitId(int set) { mGeoInitId=set; }
                /*! Returns the geo init typ */
                inline void setGeoInitType(int set) { mGeoInitType=set; }
-               /*! Returns the geo init id */
-               inline int getGeoInitId() const { return mGeoInitId; }
                /*! Returns the geo init typ */
                inline int getGeoInitType() const { return mGeoInitType; }
 
@@ -83,8 +83,8 @@ class ntlGeometryObject : public ntlGeometryClass
                inline void setGeoPartSlipValue(float set) { mGeoPartSlipValue=set; }
 
                /*! Set/get the part slip value*/
-               inline bool getOnlyThinInit() const { return mOnlyThinInit; }
-               inline void setOnlyThinInit(float set) { mOnlyThinInit=set; }
+               inline int getVolumeInit() const { return mVolumeInit; }
+               inline void setVolumeInit(int set) { mVolumeInit=set; }
 
                /*! Set/get the cast initial veocity attribute */
                void setInitialVelocity(ntlVec3Gfx set);
@@ -102,12 +102,23 @@ class ntlGeometryObject : public ntlGeometryClass
 
                /*! is the object animated? */
                inline bool getIsAnimated() const { return mIsAnimated; }
+               /*! init object anim flag */
+               bool checkIsAnimated();
+               /*! is the mesh animated? */
+               virtual bool getMeshAnimated();
+               /*! init triangle divisions */
+               virtual void calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri);
 
                /*! apply object translation at time t*/
                void applyTransformation(double t, vector<ntlVec3Gfx> *verts, vector<ntlVec3Gfx> *norms, int vstart, int vend, int forceTrafo);
 
                /*! Prepare points for moving objects */
-               void initMovingPoints(gfxReal featureSize);
+               void initMovingPoints(double time, gfxReal featureSize);
+               /*! Prepare points for animated objects */
+               void initMovingPointsAnim(
+                double srctime, vector<ntlVec3Gfx> &srcpoints,
+                double dsttime, vector<ntlVec3Gfx> &dstpoints,
+                gfxReal featureSize, ntlVec3Gfx geostart, ntlVec3Gfx geoend );
                /*! Prepare points for moving objects (copy into ret) */
                void getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3Gfx> *norms = NULL);
                /*! Calculate max. velocity on object from t1 to t2 */
@@ -126,6 +137,9 @@ class ntlGeometryObject : public ntlGeometryClass
                                vector<ntlTriangle> *triangles,
                                vector<ntlVec3Gfx>  *vertices,
                                vector<ntlVec3Gfx>  *vertNormals);
+               void sceneAddTriangleNoVert(int *trips,
+                               ntlVec3Gfx trin, bool smooth,
+                               vector<ntlTriangle> *triangles);
 
        protected:
 
@@ -144,8 +158,6 @@ class ntlGeometryObject : public ntlGeometryClass
                int mReceiveShadows;
 
                /* fluid init data */
-               /*! id of fluid init (is used in solver initialization) */
-               int mGeoInitId;
                /*! fluid object type (fluid, obstacle, accelerator etc.) */
                int mGeoInitType;
                /*! initial velocity for fluid objects */
@@ -158,7 +170,7 @@ class ntlGeometryObject : public ntlGeometryClass
                /*! part slip bc value */
                float mGeoPartSlipValue;
                /*! only init as thin object, dont fill? */
-               bool mOnlyThinInit;
+               int mVolumeInit;
 
                /*! initial offset for rot/scale */
                ntlVec3Gfx mInitialPos;
@@ -169,11 +181,13 @@ class ntlGeometryObject : public ntlGeometryClass
                
                /*! moving point/normal storage */
                vector<ntlVec3Gfx> mMovPoints;
-               vector<ntlVec3Gfx> mMovNormals;
+               // unused vector<ntlVec3Gfx> mMovNormals;
                /*! cached points for non moving objects/timeslots */
                bool mHaveCachedMov;
                vector<ntlVec3Gfx> mCachedMovPoints;
-               vector<ntlVec3Gfx> mCachedMovNormals;
+               // unused vector<ntlVec3Gfx> mCachedMovNormals;
+               /*! precomputed triangle divisions */
+               vector<int> mTriangleDivs1,mTriangleDivs2,mTriangleDivs3;
                /*! inited? */
                float mMovPntsInited;
                /*! point with max. distance from center */
index 6adc562..07be18a 100644 (file)
@@ -21,7 +21,8 @@ class ntlGeometryShader :
 
                //! Default constructor
                inline ntlGeometryShader() :
-                       ntlGeometryClass() {};
+                       ntlGeometryClass(), mOutFilename("")
+                       {};
                //! Default destructor
                virtual ~ntlGeometryShader() {};
 
@@ -42,11 +43,17 @@ class ntlGeometryShader :
                /*! notify object that dump is in progress (e.g. for field dump) */
                virtual void notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0;
 
+               /*! get ouput filename, returns global render outfile if empty */
+               std::string getOutFilename( void ) { return mOutFilename; }
+
        protected:
 
                //! vector for the objects
                vector<ntlGeometryObject *> mObjects;
 
+
+               /*! surface output name for this simulation */
+               std::string mOutFilename; 
 };
 
 #endif
index 5df9c00..c838e6f 100644 (file)
@@ -803,7 +803,7 @@ void ntlScene::buildScene(double time,bool firstInit)
 
                int vstart = mVertNormals.size();
                (*iter)->setObjectId(idCnt);
-               (*iter)->getTriangles(&mTriangles, &mVertices, &mVertNormals, idCnt);
+               (*iter)->getTriangles(time, &mTriangles, &mVertices, &mVertNormals, idCnt);
                (*iter)->applyTransformation(time, &mVertices, &mVertNormals, vstart, mVertices.size(), false );
 
                if(debugTriCollect) debMsgStd("ntlScene::buildScene",DM_MSG,"Done with "<<(*iter)->getName()<<" totTris:"<<mTriangles.size()<<" totVerts:"<<mVertices.size()<<" totNorms:"<<mVertNormals.size(), 4 );
@@ -837,7 +837,10 @@ void ntlScene::buildScene(double time,bool firstInit)
                                iter != mGeos.end(); iter++) {
                        if( (*iter)->getTypeId() & GEOCLASSTID_SHADER ) {
                                ntlGeometryShader *geoshad = (ntlGeometryShader*)(*iter);
-                               geoshad->postGeoConstrInit( mpGlob );
+                               if(geoshad->postGeoConstrInit( mpGlob )) {
+                                       errFatal("ntlScene::buildScene","Init failed for object '"<< (*iter)->getName() <<"' !", SIMWORLD_INITERROR );
+                                       return;
+                               }
                        }
                }
                mFirstInitDone = true;
index 716a620..2321cee 100644 (file)
@@ -149,8 +149,6 @@ private:
 /*! Triangle flag defines */
 #define TRI_GEOMETRY      (1<<0)
 #define TRI_CASTSHADOWS   (1<<1)
-#define TRI_MAKECAUSTICS  (1<<2)
-#define TRI_NOCAUSTICS    (1<<3)
 
 
 class ntlTriangle
@@ -217,7 +215,7 @@ private:
        //int mNormalIndex; ??
        ntlVec3Gfx mNormal;
 
-       /*! Flags for object attributes cast shadows, make caustics etc. */
+       /*! Flags for object attributes cast shadows */
        int mFlags;
 
        /*! ID of last ray that an intersection was calculated for */
index ffcaa41..2e50330 100644 (file)
@@ -201,6 +201,9 @@ private:
 
 
 
+//! global string for formatting vector output in utilities.cpp
+extern char *globVecFormatStr;
+
 /*************************************************************************
   Outputs the object in human readable form using the format
   [x,y,z]
@@ -209,7 +212,10 @@ template<class Scalar>
 std::ostream&
 operator<<( std::ostream& os, const ntlVector3Dim<Scalar>& i )
 {
-  os << '[' << i[0] << ", " << i[1] << ", " << i[2] << ']';
+       char buf[256];
+       snprintf(buf,256,globVecFormatStr,i[0],i[1],i[2]);
+       os << std::string(buf); 
+  //os << '[' << i[0] << ", " << i[1] << ", " << i[2] << ']';
   return os;
 }
 
index 8d1d5d4..6b4dc7b 100644 (file)
@@ -34,6 +34,11 @@ void setPointers( ntlRenderGlobals *setglob);
 /******************************************************************************
  * Constructor
  *****************************************************************************/
+
+ntlWorld::ntlWorld() {
+       initDefaults();
+}
+
 ntlWorld::ntlWorld(string filename, bool commandlineMode) 
 {
 #ifndef ELBEEM_PLUGIN
@@ -60,12 +65,18 @@ ntlWorld::ntlWorld(string filename, bool commandlineMode)
 #endif // ELBEEM_PLUGIN
 }
 
-ntlWorld::ntlWorld(elbeemSimulationSettings *settings)
+
+int globalDomainCounter = 1;
+int ntlWorld::addDomain(elbeemSimulationSettings *settings)
 {
-       initDefaults();
-       // todo init settings
+       // create domain obj
        SimulationObject *sim = new SimulationObject();
+       char simname[100];
+       snprintf(simname,100,"domain%04d",globalDomainCounter);
+       globalDomainCounter++;
+       sim->setName(std::string(simname));
        mpGlob->getSims()->push_back( sim );
+
        // important - add to both, only render scene objects are free'd 
        mpGlob->getRenderScene()->addGeoClass( sim );
        mpGlob->getSimScene()->addGeoClass( sim );
@@ -80,10 +91,7 @@ ntlWorld::ntlWorld(elbeemSimulationSettings *settings)
        Parametrizer *param = sim->getParametrizer();
        param->setSize( settings->resolutionxyz );
        param->setDomainSize( settings->realsize );
-       param->setViscosity( settings->viscosity );
-       param->setGravity( ParamVec(settings->gravity[0], settings->gravity[1], settings->gravity[2]) );
        param->setAniStart( settings->animStart );
-       param->setAniFrameTimeChannel( settings->aniFrameTime );
        param->setNormalizedGStar( settings->gstar );
 
        // init domain channels
@@ -98,21 +106,30 @@ ntlWorld::ntlWorld(elbeemSimulationSettings *settings)
        valv.clear(); time.clear(); elbeemSimplifyChannelVec3(channel,&size); \
        for(int i=0; i<size; i++) { valv.push_back( ParamVec(channel[4*i+0],channel[4*i+1],channel[4*i+2]) ); time.push_back( channel[4*i+3] ); } 
 
-       INIT_CHANNEL_FLOAT(settings->channelViscosity, settings->channelSizeViscosity);
-       param->initViscosityChannel(valf,time);
+       param->setViscosity( settings->viscosity );
+       if((settings->channelViscosity)&&(settings->channelSizeViscosity>0)) {
+               INIT_CHANNEL_FLOAT(settings->channelViscosity, settings->channelSizeViscosity);
+               param->initViscosityChannel(valf,time); }
 
-       INIT_CHANNEL_VEC(settings->channelGravity, settings->channelSizeGravity);
-       param->initGravityChannel(valv,time);
+       param->setGravity( ParamVec(settings->gravity[0], settings->gravity[1], settings->gravity[2]) );
+       if((settings->channelGravity)&&(settings->channelSizeGravity>0)) {
+               INIT_CHANNEL_VEC(settings->channelGravity, settings->channelSizeGravity);
+               param->initGravityChannel(valv,time); }
 
-       INIT_CHANNEL_FLOAT(settings->channelFrameTime, settings->channelSizeFrameTime);
-       param->initAniFrameTimeChannel(valf,time);
+       param->setAniFrameTimeChannel( settings->aniFrameTime );
+       if((settings->channelFrameTime)&&(settings->channelSizeFrameTime>0)) {
+               INIT_CHANNEL_FLOAT(settings->channelFrameTime, settings->channelSizeFrameTime);
+               param->initAniFrameTimeChannel(valf,time); }
 
 #undef INIT_CHANNEL_FLOAT
 #undef INIT_CHANNEL_VEC
        
-       mpGlob->setAniFrames( settings->noOfFrames );
+       // might be set by previous domain
+       if(mpGlob->getAniFrames() < settings->noOfFrames)       mpGlob->setAniFrames( settings->noOfFrames );
+       // set additionally to SimulationObject->mOutFilename
        mpGlob->setOutFilename( settings->outputPath );
-       // further init in postGeoConstrInit/initializeLbmSimulation of SimulationObject
+
+       return 0;
 }
 
 void ntlWorld::initDefaults()
@@ -149,18 +166,19 @@ void ntlWorld::initDefaults()
 
 void ntlWorld::finishWorldInit()
 {
-       if(!SIMWORLD_OK()) return;
+       if(! isSimworldOk() ) return;
 
        // init the scene for the first time
   long sstartTime = getTime();
 
        // first init sim scene for geo setup
        mpGlob->getSimScene()->buildScene(0.0, true);
+       if(! isSimworldOk() ) return;
        mpGlob->getRenderScene()->buildScene(0.0, true);
+       if(! isSimworldOk() ) return;
        long sstopTime = getTime();
        debMsgStd("ntlWorld::ntlWorld",DM_MSG,"Scene build time: "<< getTimeString(sstopTime-sstartTime) <<" ", 10);
 
-       if(!SIMWORLD_OK()) return;
        // TODO check simulations, run first steps
        mFirstSim = -1;
        if(mpSims->size() > 0) {
@@ -204,6 +222,9 @@ void ntlWorld::finishWorldInit()
                        if(!mpGlob->getSingleFrameMode()) debMsgStd("ntlWorld::ntlWorld",DM_WARNING,"No active simulations!", 1);
                }
        }
+
+       if(! isSimworldOk() ) return;
+       setElbeemState( SIMWORLD_INITED );
 }
 
 
@@ -237,11 +258,9 @@ void ntlWorld::setSingleFrameOut(string singleframeFilename) {
  *****************************************************************************/
 
 // blender interface
-#if ELBEEM_BLENDER==1
-extern "C" {
-       void simulateThreadIncreaseFrame(void);
-}
-#endif // ELBEEM_BLENDER==1
+//#if ELBEEM_BLENDER==1
+//extern "C" { void simulateThreadIncreaseFrame(void); }
+//#endif // ELBEEM_BLENDER==1
 
 int ntlWorld::renderAnimation( void )
 {
@@ -258,20 +277,32 @@ int ntlWorld::renderAnimation( void )
        } 
 
        mThreadRunning = true; // not threaded, but still use the same flags
-       renderScene();
+       if(getElbeemState() == SIMWORLD_INITED) {
+               renderScene();
+       } else if(getElbeemState() == SIMWORLD_STOP) {
+               // dont render now, just continue
+               setElbeemState( SIMWORLD_INITED );
+               mFrameCnt--; // counted one too many from last abort...
+       } else {
+               debMsgStd("ntlWorld::renderAnimation",DM_NOTIFY,"Not properly inited, stopping...",1);
+               return 1;
+       }
+       
        if(mpSims->size() <= 0) {
                debMsgStd("ntlWorld::renderAnimation",DM_NOTIFY,"No simulations found, stopping...",1);
                return 1;
        }
 
-       for(mFrameCnt=0; ((mFrameCnt<mpGlob->getAniFrames()) && (!getStopRenderVisualization() )); mFrameCnt++) {
+       bool simok = true;
+       for( ; ((mFrameCnt<mpGlob->getAniFrames()) && (!getStopRenderVisualization() ) && (simok)); mFrameCnt++) {
                if(!advanceSims(mFrameCnt)) {
                        renderScene();
-#if ELBEEM_BLENDER==1
+//#if ELBEEM_BLENDER==1
                        // update Blender gui display after each frame
-                       simulateThreadIncreaseFrame();
-#endif // ELBEEM_BLENDER==1
+                       //simulateThreadIncreaseFrame();
+//#endif // ELBEEM_BLENDER==1
                } // else means sim panicked, so dont render...
+               else { simok=false; }
        }
        mThreadRunning = false;
        return 0;
@@ -285,8 +316,10 @@ int ntlWorld::renderAnimation( void )
 int ntlWorld::renderVisualization( bool multiThreaded ) 
 {
 #ifndef NOGUI
-       //gfxReal deltat = 0.0015;
+       if(getElbeemState() != SIMWORLD_INITED) { return 0; }
+
        if(multiThreaded) mThreadRunning = true;
+       // TODO, check global state?
        while(!getStopRenderVisualization()) {
 
                if(mpSims->size() <= 0) {
@@ -315,7 +348,7 @@ int ntlWorld::renderVisualization( bool multiThreaded )
                                warnMsg("ntlWorld::advanceSims","All sims panicked... stopping thread" );
                                setStopRenderVisualization( true );
                        }
-                       if(!SIMWORLD_OK()) {
+                       if(! isSimworldOk() ) {
                                warnMsg("ntlWorld::advanceSims","World state error... stopping" );
                                setStopRenderVisualization( true );
                        }
@@ -368,7 +401,11 @@ int ntlWorld::advanceSims(int framenum)
 {
        bool done = false;
        bool allPanic = true;
-       //debMsgStd("ntlWorld::advanceSims",DM_MSG,"Advancing sims to "<<targetTime, 10 ); // timedebug
+
+       // stop/quit, dont display/render
+       if(getElbeemState()==SIMWORLD_STOP) { 
+               return 1;
+       }
 
        for(size_t i=0;i<mpSims->size();i++) { (*mpSims)[i]->setFrameNum(framenum); }
        double targetTime = mSimulationTime + (*mpSims)[mFirstSim]->getFrameTime(framenum);
@@ -378,16 +415,15 @@ int ntlWorld::advanceSims(int framenum)
                done=true; allPanic=false; 
        }
 
-#if ELBEEM_BLENDER==1
+       int gstate = 0;
+//#if ELBEEM_BLENDER==1
        // same as solver_main check, but no mutex check here
-       if(getGlobalBakeState()<0) {
-               // this means abort... cause panic
-               allPanic = true; done = true;
-       }
-#endif // ELBEEM_BLENDER==1
+       //gstate = getGlobalBakeState();
+       //if(gstate<0) { allPanic = true; done = true; } // this means abort... cause panic
+//#endif // ELBEEM_BLENDER==1
 
        // step all the sims, and check for panic
-       debMsgStd("ntlWorld::advanceSims",DM_MSG, " sims "<<mpSims->size()<<" t"<<targetTime<<" done:"<<done<<" panic:"<<allPanic, 10); // debug // timedebug
+       debMsgStd("ntlWorld::advanceSims",DM_MSG, " sims "<<mpSims->size()<<" t"<<targetTime<<" done:"<<done<<" panic:"<<allPanic<<" gstate:"<<gstate, 10); // debug // timedebug
        while(!done) {
                double nextTargetTime = (*mpSims)[mFirstSim]->getCurrentTime() + (*mpSims)[mFirstSim]->getTimestep();
                singleStepSims(nextTargetTime);
@@ -395,7 +431,7 @@ int ntlWorld::advanceSims(int framenum)
                // check target times
                done = true;
                allPanic = false;
-               //if((framenum>0) && (nextTargetTime<=(*mpSims)[mFirstSim]->getCurrentTime()) ) { 
+               
                if((*mpSims)[mFirstSim]->getTimestep() <1e-9 ) { 
                        // safety check, avoid timesteps that are too small
                        errMsg("ntlWorld::advanceSims","Invalid time step, causing panic! curr:"<<(*mpSims)[mFirstSim]->getCurrentTime()<<" next:"<<nextTargetTime<<", stept:"<< (*mpSims)[mFirstSim]->getTimestep() );
@@ -468,12 +504,11 @@ void ntlWorld::singleStepSims(double targetTime) {
 int ntlWorld::renderScene( void )
 {
 #ifndef ELBEEM_PLUGIN
-       char nrStr[5];                                                                                                          /* nr conversion */
-       //std::ostringstream outfilename("");                                     /* ppm file */
-       std::ostringstream outfn_conv("");                                              /* converted ppm with other suffix */
-  ntlRenderGlobals *glob;                      /* storage for global rendering parameters */
-  myTime_t timeStart,totalStart,timeEnd;               /* measure user running time */
-  myTime_t rendStart,rendEnd;                          /* measure user rendering time */
+       char nrStr[5];                                                                                                          // nr conversion 
+       std::ostringstream outfn_conv("");                      // converted ppm with other suffix 
+  ntlRenderGlobals *glob;                      // storage for global rendering parameters 
+  myTime_t timeStart,totalStart,timeEnd;               // measure user running time 
+  myTime_t rendStart,rendEnd;                          // measure user rendering time 
   glob = mpGlob;
 
        /* check if picture already exists... */
index 131f182..18f388a 100644 (file)
@@ -23,16 +23,18 @@ class ntlRandomStream;
 class ntlWorld
 {
        public:
+               /*! Constructor for API init */
+               ntlWorld();
                /*! Constructor */
                ntlWorld(string filename, bool commandlineMode);
-               /*! Constructor for API init */
-               ntlWorld(elbeemSimulationSettings *simSettings);
                /*! Destructor */
                virtual ~ntlWorld( void );
                /*! default init for all contructors */
                void initDefaults();
                /*! common world contruction stuff once the scene is set up */
                void finishWorldInit();
+               /*! add domain for API init */
+               int addDomain(elbeemSimulationSettings *simSettings);
 
                /*! render a whole animation (command line mode) */
                int renderAnimation( void );
index 50d544d..b9dc1c3 100644 (file)
@@ -55,7 +55,6 @@ Parametrizer::Parametrizer( void ) :
        mSizex(50), mSizey(50), mSizez(50),
        mTimeFactor( 1.0 ),
        mcAniFrameTime(0.0001),
-       mAniFrameTime(0.0001), 
        mTimeStepScale(1.0),
        mAniStart(0.0),
        mExtent(1.0, 1.0, 1.0), //mSurfaceTension( 0.0 ),
@@ -149,17 +148,11 @@ void Parametrizer::parseAttrList()
  *! advance to next render/output frame 
  *****************************************************************************/
 void Parametrizer::setFrameNum(int frame) {
-       //double oldval = mAniFrameTime;
-       //mAniFrameTime = mcAniFrameTime.get(frametime);
-       //if(mAniFrameTime<0.0) {
-               //errMsg("Parametrizer::setFrameNum","Invalid frame time:"<<mAniFrameTime<<" at frame "<<frame<<", resetting to "<<oldval);
-               //mAniFrameTime = oldval; }
-       //errMsg("ChannelAnimDebug","anim: anif"<<mAniFrameTime<<" at "<<frame<<" ");
-       // debug getAttributeList()->find("p_aniframetime")->print();
        mFrameNum = frame;
        if(DEBUG_PARAMCHANNELS) errMsg("DEBUG_PARAMCHANNELS","setFrameNum frame-num="<<mFrameNum);
 }
 /*! get time of an animation frame (renderer)  */
+// also used by: mpParam->getCurrentAniFrameTime() , e.g. for velocity dump
 ParamFloat Parametrizer::getAniFrameTime( int frame )   { 
        double frametime = (double)frame;
        ParamFloat anift = mcAniFrameTime.get(frametime);
@@ -168,7 +161,7 @@ ParamFloat Parametrizer::getAniFrameTime( int frame )   {
                errMsg("Parametrizer::setFrameNum","Invalid frame time:"<<anift<<" at frame "<<frame<<", resetting to "<<resetv);
                anift = resetv; 
        }
-       if(DEBUG_PARAMCHANNELS) errMsg("DEBUG_PARAMCHANNELS","getAniFrameTime frame="<<frame<<", frametime="<<anift<<" ");
+       if((0)|| (DEBUG_PARAMCHANNELS)) errMsg("DEBUG_PARAMCHANNELS","getAniFrameTime frame="<<frame<<", frametime="<<anift<<" ");
        return anift; 
 }
 
@@ -247,7 +240,7 @@ int Parametrizer::calculateAniStepsPerFrame(int frame)   {
        }
        int value = (int)(getAniFrameTime(frame)/mTimestep); 
        if((value<0) || (value>1000000)) {
-               errFatal("Parametrizer::calculateAniStepsPerFrame", "Invalid step-time (="<<mAniFrameTime<<") <> ani-frame-time ("<<mTimestep<<") settings, aborting...", SIMWORLD_INITERROR);
+               errFatal("Parametrizer::calculateAniStepsPerFrame", "Invalid step-time (="<<value<<") <> ani-frame-time ("<<mTimestep<<") settings, aborting...", SIMWORLD_INITERROR);
                return 1;
        }
        return value;
@@ -373,7 +366,7 @@ errMsg("Warning","Used z-dir for gstar!");
                        maxDeltaT = sqrt( gStar*mCellSize *mTimeStepScale /forceStrength );
                } else {
                        // use 1 lbm setp = 1 anim step as max
-                       maxDeltaT = mAniFrameTime;
+                       maxDeltaT = getAniFrameTime(0);
                }
 
                if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," targeted step time = "<<maxDeltaT, 10);
@@ -436,16 +429,16 @@ errMsg("Warning","Used z-dir for gstar!");
                
                if(!checkSeenValues(PARAM_ANIFRAMETIME)) {
                        errFatal("Parametrizer::calculateAllMissingValues"," Warning no ani frame time given!", SIMWORLD_INITERROR);
-                       mAniFrameTime = mTimestep;
+                       setAniFrameTimeChannel( mTimestep );
                } 
-               //mAniFrameTime = mAniFrames * mTimestep;
+
                if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," ani frame steps = "<<calculateAniStepsPerFrame(mFrameNum)<<" for frame "<<mFrameNum, 1);
 
                if((checkSeenValues(PARAM_ANISTART))&&(calculateAniStart()>0)) {
                        if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," ani start steps = "<<calculateAniStart()<<" ",1); 
                }
                        
-               if(!SIMWORLD_OK()) return false;
+               if(! isSimworldOk() ) return false;
                // everything ok
                return true;
        }
@@ -471,6 +464,51 @@ errMsg("Warning","Used z-dir for gstar!");
 }
 
 
+/******************************************************************************
+ * init debug functions
+ *****************************************************************************/
+
+
+/*! set kinematic viscosity */
+void Parametrizer::setViscosity(ParamFloat set) { 
+       mcViscosity = AnimChannel<ParamFloat>(set); 
+       seenThis( PARAM_VISCOSITY ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcViscosity set = "<< mcViscosity.printChannel() ); }
+}
+void Parametrizer::initViscosityChannel(vector<ParamFloat> val, vector<double> time) { 
+       mcViscosity = AnimChannel<ParamFloat>(val,time); 
+       seenThis( PARAM_VISCOSITY ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcViscosity initc = "<< mcViscosity.printChannel() ); }
+}
+
+/*! set the external force */
+void Parametrizer::setGravity(ParamFloat setx, ParamFloat sety, ParamFloat setz) { 
+       mcGravity = AnimChannel<ParamVec>(ParamVec(setx,sety,setz)); 
+       seenThis( PARAM_GRAVITY ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcGravity set = "<< mcGravity.printChannel() ); }
+}
+void Parametrizer::setGravity(ParamVec set) { 
+       mcGravity = AnimChannel<ParamVec>(set); 
+       seenThis( PARAM_GRAVITY ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcGravity set = "<< mcGravity.printChannel() ); }
+}
+void Parametrizer::initGravityChannel(vector<ParamVec> val, vector<double> time) { 
+       mcGravity = AnimChannel<ParamVec>(val,time); 
+       seenThis( PARAM_GRAVITY ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcGravity initc = "<< mcGravity.printChannel() ); }
+}
+
+/*! set time of an animation frame (renderer)  */
+void Parametrizer::setAniFrameTimeChannel(ParamFloat set) { 
+       mcAniFrameTime = AnimChannel<ParamFloat>(set); 
+       seenThis( PARAM_ANIFRAMETIME ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcAniFrameTime set = "<< mcAniFrameTime.printChannel() ); }
+}
+void Parametrizer::initAniFrameTimeChannel(vector<ParamFloat> val, vector<double> time) { 
+       mcAniFrameTime = AnimChannel<ParamFloat>(val,time); 
+       seenThis( PARAM_ANIFRAMETIME ); 
+       if(DEBUG_PARAMCHANNELS) { errMsg("DebugChannels","Parametrizer::mcAniFrameTime initc = "<< mcAniFrameTime.printChannel() ); }
+}
 
 // OLD interface stuff
 // reactivate at some point?
index 471fb57..c0df399 100644 (file)
@@ -119,22 +119,25 @@ class Parametrizer {
                /*! calculate real world [m/s] velocity from lattice value */
                ParamVec calculateRwVelocityFromLatt( ParamVec ivel );
 
-
-               /*! set kinematic viscosity */
-               void setViscosity(ParamFloat set) { mcViscosity = AnimChannel<ParamFloat>(set); seenThis( PARAM_VISCOSITY ); }
-               void initViscosityChannel(vector<ParamFloat> val, vector<double> time) { mcViscosity = AnimChannel<ParamFloat>(val,time); }
-
                /*! set speed of sound */
                void setSoundSpeed(ParamFloat set) { mSoundSpeed = set; seenThis( PARAM_SOUNDSPEED ); }
                /*! get speed of sound */
                ParamFloat getSoundSpeed( void )   { return mSoundSpeed; }
 
+               /*! set kinematic viscosity */
+               void setViscosity(ParamFloat set);
+               void initViscosityChannel(vector<ParamFloat> val, vector<double> time);
+
                /*! set the external force */
-               void setGravity(ParamFloat setx, ParamFloat sety, ParamFloat setz) { mcGravity = AnimChannel<ParamVec>(ParamVec(setx,sety,setz)); seenThis( PARAM_GRAVITY ); }
-               void setGravity(ParamVec set) { mcGravity = AnimChannel<ParamVec>(set); seenThis( PARAM_GRAVITY ); }
-               void initGravityChannel(vector<ParamVec> val, vector<double> time) { mcGravity = AnimChannel<ParamVec>(val,time); }
+               void setGravity(ParamFloat setx, ParamFloat sety, ParamFloat setz);
+               void setGravity(ParamVec set);
+               void initGravityChannel(vector<ParamVec> val, vector<double> time);
                ParamVec getGravity(double time) { return mcGravity.get( time ); }
 
+               /*! set time of an animation frame (renderer)  */
+               void setAniFrameTimeChannel(ParamFloat set);
+               void initAniFrameTimeChannel(vector<ParamFloat> val, vector<double> time);
+
                /*! set the length of a single time step */
                void setTimestep(ParamFloat set) { mTimestep = set; seenThis( PARAM_STEPTIME ); }
                /*! get the length of a single time step */
@@ -155,11 +158,6 @@ class Parametrizer {
                void setSize(int ijk)            { mSizex = ijk; mSizey = ijk; mSizez = ijk; seenThis( PARAM_SIZE ); }
                void setSize(int i,int j, int k) { mSizex = i; mSizey = j; mSizez = k; seenThis( PARAM_SIZE ); }
 
-               /*! set time of an animation frame (renderer)  */
-               //void setAniFrameTime(ParamFloat set) { mAniFrameTime = set; seenThis( PARAM_ANIFRAMETIME ); }
-               void setAniFrameTimeChannel(ParamFloat set) { mcAniFrameTime = AnimChannel<ParamFloat>(set); seenThis( PARAM_ANIFRAMETIME ); }
-               void initAniFrameTimeChannel(vector<ParamFloat> val, vector<double> time) { mcAniFrameTime = AnimChannel<ParamFloat>(val,time); seenThis( PARAM_ANIFRAMETIME ); }
-
                /*! set starting time of the animation (renderer) */
                void setAniStart(ParamFloat set) { mAniStart = set; seenThis( PARAM_ANISTART ); }
                /*! get starting time of the animation (renderer) */
@@ -266,8 +264,6 @@ class Parametrizer {
 
                /*! for renderer - length of an animation step [s] */
                AnimChannel<ParamFloat> mcAniFrameTime;
-               /*! current value for next frame */
-               ParamFloat mAniFrameTime;
                /*! time step scaling factor for testing/debugging */
                ParamFloat mTimeStepScale;
 
index 8286746..c49c1b1 100644 (file)
@@ -33,23 +33,17 @@ ParticleTracer::ParticleTracer() :
        mPartSize(0.01),
        mStart(-1.0), mEnd(1.0),
        mSimStart(-1.0), mSimEnd(1.0),
-       mPartScale(1.0) , mPartHeadDist( 0.5 ), mPartTailDist( -4.5 ), mPartSegments( 4 ),
+       mPartScale(0.1) , mPartHeadDist( 0.1 ), mPartTailDist( -0.1 ), mPartSegments( 4 ),
        mValueScale(0),
        mValueCutoffTop(0.0), mValueCutoffBottom(0.0),
-       mDumpParts(0), mDumpText(0), mDumpTextFile(""), mShowOnly(0), 
-       mNumInitialParts(0), mpTrafo(NULL)
+       mDumpParts(0), //mDumpText(0), 
+       mDumpTextFile(""), 
+       mDumpTextInterval(0.), mDumpTextLastTime(0.), mDumpTextCount(0),
+       mShowOnly(0), 
+       mNumInitialParts(0), mpTrafo(NULL),
+       mInitStart(-1.), mInitEnd(-1.),
+       mPrevs(), mTrailTimeLast(0.), mTrailInterval(-1.), mTrailLength(0)
 {
-       // check env var
-//#ifdef ELBEEM_PLUGIN
-       //mDumpParts=1; // default on
-//#endif // ELBEEM_PLUGIN
-       //if(getenv("ELBEEM_DUMPPARTICLE")) { // DEBUG!
-               //int set = atoi( getenv("ELBEEM_DUMPPARTICLE") );
-               //if((set>=0)&&(set!=mDumpParts)) {
-                       //mDumpParts=set;
-                       //debMsgStd("ParticleTrace",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mDumpParts to "<<set<<","<<mDumpParts,8);
-               //}
-       //}
        debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10);
 };
 
@@ -66,7 +60,7 @@ void ParticleTracer::parseAttrList(AttributeList *att)
        mpAttrs = att;
 
        mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false);
-       errMsg(" NUMP"," "<<mNumInitialParts);
+       //errMsg(" NUMP"," "<<mNumInitialParts);
        mPartScale    = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false);
        mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false);
        mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false);
@@ -76,19 +70,23 @@ void ParticleTracer::parseAttrList(AttributeList *att)
        mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false);
 
        mDumpParts   = mpAttrs->readInt  ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false);
-       mDumpText    = mpAttrs->readInt  ("part_textdump",mDumpText, "ParticleTracer","mDumpText", false);
+       // mDumpText deprecatd, use mDumpTextInterval>0. instead
        mShowOnly    = mpAttrs->readInt  ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false);
        mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false);
+       mDumpTextInterval= mpAttrs->readFloat("part_textdumpinterval",mDumpTextInterval, "ParticleTracer","mDumpTextInterval", false);
 
        string matPart;
        matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false);
        setMaterialName( matPart );
 
+       mInitStart = mpAttrs->readFloat("part_initstart",mInitStart, "ParticleTracer","mInitStart", false);
+       mInitEnd   = mpAttrs->readFloat("part_initend",  mInitEnd, "ParticleTracer","mInitEnd", false);
+
        // unused...
-       int mTrailLength  = 0; // UNUSED
-       int mTrailInterval= 0; // UNUSED
+       //int mTrailLength  = 0; // UNUSED
+       //int mTrailInterval= 0; // UNUSED
        mTrailLength  = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false);
-       mTrailInterval= mpAttrs->readInt("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
+       mTrailInterval= mpAttrs->readFloat("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false);
 
        // restore old list
        mpAttrs = tempAtt;
@@ -134,24 +132,14 @@ void ParticleTracer::addParticle(float x, float y, float z)
 {
        ntlVec3Gfx p(x,y,z);
        ParticleObject part( p );
-       //mParts.push_back( part );
-       // TODO handle other arrays?
-       //part.setActive( false );
        mParts.push_back( part );
-       //for(size_t l=0; l<mParts.size(); l++) {
-               // add deactivated particles to other arrays
-               // deactivate further particles
-               //if(l>1) {
-                       //mParts[l][ mParts.size()-1 ].setActive( false );
-               //}
-       //}
 }
 
 
 void ParticleTracer::cleanup() {
        // cleanup
        int last = (int)mParts.size()-1;
-       if(mDumpText>0) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }
+       if(mDumpTextInterval>0.) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; }
 
        for(int i=0; i<=last; i++) {
                if( mParts[i].getActive()==false ) {
@@ -166,7 +154,7 @@ void ParticleTracer::cleanup() {
  *! dump particles if desired 
  *****************************************************************************/
 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
-       debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr, 10); // DEBUG
+       debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts, 10); // DEBUG
 
        if(
                        (dumptype==DUMP_FULLGEOMETRY)&&
@@ -186,13 +174,11 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
                        numParts = 0;
                        for(size_t i=0; i<mParts.size(); i++) {
                                if(!mParts[i].getActive()) continue;
-                               //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
                                numParts++;
                        }
                        gzwrite(gzf, &numParts, sizeof(numParts));
                        for(size_t i=0; i<mParts.size(); i++) {
                                if(!mParts[i].getActive()) { continue; }
-                               //if(mParts[i].getLifeTime()<30) { continue; } //? CHECK
                                ParticleObject *p = &mParts[i];
                                //int type = p->getType();  // export whole type info
                                int type = p->getFlags(); // debug export whole type & status info
@@ -221,17 +207,21 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
                        gzclose( gzf );
                }
        } // dump?
+}
 
+void ParticleTracer::checkDumpTextPositions(double simtime) {
        // dfor partial & full dump
-       if(mDumpText>0) {
+       errMsg("ParticleTracer::checkDumpTextPositions","t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval);
+
+       if((mDumpTextInterval>0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) {
                // dump to binary file
                std::ostringstream boutfilename("");
                if(mDumpTextFile.length()>1) {   
                        boutfilename << mDumpTextFile <<  ".cpart2"; 
                } else {                           
-                       boutfilename << outfilename <<"_particles" <<  ".cpart2"; 
+                       boutfilename << boutfilename <<"_particles" <<  ".cpart2"; 
                }
-               debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7);
+               debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" " , 7);
 
                int numParts = 0;
                // only dump bubble particles
@@ -242,27 +232,34 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
                }
 
                // output to text file
-               gzFile gzf;
-               if(frameNr==0) {
-                       gzf = gzopen(boutfilename.str().c_str(), "w0");
+               //gzFile gzf;
+               FILE *stf;
+               if(mDumpTextCount==0) {
+                       //gzf = gzopen(boutfilename.str().c_str(), "w0");
+                       stf = fopen(boutfilename.str().c_str(), "w");
 
-                       gzprintf( gzf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
+                       fprintf( stf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts);
                        // fixed time scale for now
-                       gzprintf( gzf, "T %f \n\n", 1.0);
+                       fprintf( stf, "T %f \n\n", 1.0);
                } else {
-                       gzf = gzopen(boutfilename.str().c_str(), "a+0");
+                       //gzf = gzopen(boutfilename.str().c_str(), "a+0");
+                       stf = fopen(boutfilename.str().c_str(), "a+");
                }
 
                // add current set
-               if(gzf) {
-                       gzprintf( gzf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", frameNr, simtime, numParts );
-                       gzprintf( gzf, "S %f \n\n", simtime );
+               if(stf) {
+                       fprintf( stf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", mDumpTextCount, simtime, numParts );
+                       fprintf( stf, "S %f \n\n", simtime );
                        
                        for(size_t i=0; i<mParts.size(); i++) {
                                ParticleObject *p = &mParts[i];
                                ntlVec3Gfx pos = p->getPos();
                                float size = p->getSize();
-                               if(!mParts[i].getActive()) { size=0.; } // switch "off"
+                               float infl = 1.;
+                               //if(!mParts[i].getActive()) { size=0.; } // switch "off"
+                               if(!mParts[i].getActive()) { infl=0.; } // switch "off"
+                               if(!mParts[i].getInFluid()) { infl=0.; } // switch "off"
+                               if(mParts[i].getLifeTime()<0.) { infl=0.; } // not yet active...
 
                                pos = (*mpTrafo) * pos;
                                ntlVec3Gfx v = p->getVel();
@@ -270,23 +267,45 @@ void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,str
                                v[1] *= mpTrafo->value[1][1];
                                v[2] *= mpTrafo->value[2][2];
                                
-                               gzprintf( gzf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
-                               gzprintf( gzf, "s %f \n", size );
-                               gzprintf( gzf, "\n", size );
+                               fprintf( stf, "P %f %f %f \n", pos[0],pos[1],pos[2] );
+                               if(size!=1.0) fprintf( stf, "s %f \n", size );
+                               if(infl!=1.0) fprintf( stf, "i %f \n", infl );
+                               fprintf( stf, "\n" );
                        }
 
-                       gzprintf( gzf, "# %d end  ", frameNr );
-                       gzclose( gzf );
+                       fprintf( stf, "# %d end  ", mDumpTextCount );
+                       //gzclose( gzf );
+                       fclose( stf );
+
+                       mDumpTextCount++;
                }
+
+               mDumpTextLastTime += mDumpTextInterval;
        }
+
 }
 
 
+void ParticleTracer::checkTrails(double time) {
+       if(mTrailLength<1) return;
+       if(time-mTrailTimeLast > mTrailInterval) {
+
+               if( (int)mPrevs.size() < mTrailLength) mPrevs.resize( mTrailLength );
+               for(int i=mPrevs.size()-1; i>0; i--) {
+                       mPrevs[i] = mPrevs[i-1];
+                       //errMsg("TRAIL"," from "<<i<<" to "<<(i-1) );
+               }
+               mPrevs[0] = mParts;
+
+               mTrailTimeLast += mTrailInterval;
+       }
+}
+
 
 /******************************************************************************
  * Get triangles for rendering
  *****************************************************************************/
-void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles, 
+void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles, 
                                                                                                         vector<ntlVec3Gfx> *vertices, 
                                                                                                         vector<ntlVec3Gfx> *normals, int objectId )
 {
@@ -295,6 +314,7 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
        vertices = NULL; triangles = NULL;
        normals = NULL; objectId = 0;
 #else // ELBEEM_PLUGIN
+       int pcnt = 0;
        // currently not used in blender
        objectId = 0; // remove, deprecated
        if(mDumpParts>1) { 
@@ -303,21 +323,28 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
 
        const bool debugParts = false;
        int tris = 0;
-       gfxReal partNormSize = 0.01 * mPartScale;
        int segments = mPartSegments;
-
-       //int lnewst = mTrailLength-1;
-       // TODO get rid of mPart[X] array
-       //? int loldst = mTrailLength-2;
-       // trails gehen nicht so richtig mit der
-       // richtung der partikel...
-
        ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2]));
        ntlVec3Gfx trans = mStart;
+       t = 0.; // doesnt matter
 
-       for(size_t i=0; i<mParts.size(); i++) {
-               //mParts[0][i].setActive(true);
-               ParticleObject *p = &mParts[i];
+       for(size_t t=0; t<mPrevs.size()+1; t++) {
+               vector<ParticleObject> *dparts;
+               if(t==0) {
+                       dparts = &mParts;
+               } else {
+                       dparts = &mPrevs[t-1];
+               }
+               //errMsg("TRAILT","prevs"<<t<<"/"<<mPrevs.size()<<" parts:"<<dparts->size() );
+
+       gfxReal partscale = mPartScale;
+       if(t>1) { 
+               partscale *= (gfxReal)(mPrevs.size()+1-t) / (gfxReal)(mPrevs.size()+1); 
+       }
+       gfxReal partNormSize = 0.01 * partscale;
+       //for(size_t i=0; i<mParts.size(); i++) {
+       for(size_t i=0; i<dparts->size(); i++) {
+               ParticleObject *p = &( (*dparts)[i] ); //  mParts[i];
 
                if(mShowOnly!=10) {
                        // 10=show only deleted
@@ -339,34 +366,36 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
                        if(type&PART_INTER) continue;
                }
 
+               pcnt++;
                ntlVec3Gfx pnew = p->getPos();
                if(type&PART_FLOAT) { // WARNING same handling for dump!
+                       if(p->getStatus()&PART_IN) { pnew[2] += 0.8; } // offset for display
                        // add one gridcell offset
                        //pnew[2] += 1.0; 
                }
+#if LBMDIM==2
+               pnew[2] += 0.001; // DEBUG
+               pnew[2] += 0.009; // DEBUG
+#endif 
 
                ntlVec3Gfx pdir = p->getVel();
                gfxReal plen = normalize( pdir );
                if( plen < 1e-05) pdir = ntlVec3Gfx(-1.0 ,0.0 ,0.0);
                ntlVec3Gfx pos = (*mpTrafo) * pnew;
-               //ntlVec3Gfx pos = pnew; // T
-               //pos[0] = pos[0]*scale[0]+trans[0]; // T
-               //pos[1] = pos[1]*scale[1]+trans[1]; // T
-               //pos[2] = pos[2]*scale[2]+trans[2]; // T
                gfxReal partsize = 0.0;
                if(debugParts) errMsg("DebugParts"," i"<<i<<" new"<<pnew<<" vel"<<pdir<<"   pos="<<pos );
                //if(i==0 &&(debugParts)) errMsg("DebugParts"," i"<<i<<" new"<<pnew[0]<<" pos="<<pos[0]<<" scale="<<scale[0]<<"  t="<<trans[0] );
                
                // value length scaling?
                if(mValueScale==1) {
-                       partsize = mPartScale * plen;
+                       partsize = partscale * plen;
                } else if(mValueScale==2) {
                        // cut off scaling
                        if(plen > mValueCutoffTop) continue;
                        if(plen < mValueCutoffBottom) continue;
-                       partsize = mPartScale * plen;
+                       partsize = partscale * plen;
                } else {
-                       partsize = mPartScale; // no length scaling
+                       partsize = partscale; // no length scaling
                }
                //if(type&(PART_DROP|PART_BUBBLE)) 
                partsize *= p->getSize()/5.0;
@@ -382,11 +411,7 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
                ntlVec3Gfx cv1 = pdir;
                ntlVec3Gfx cv2 = ntlVec3Gfx(pdir[1], -pdir[0], 0.0);
                ntlVec3Gfx cv3 = cross( cv1, cv2);
-               for(int l=0; l<3; l++) {
-                       //? cvmat.value[l][0] = cv1[l];
-                       //? cvmat.value[l][1] = cv2[l];
-                       //? cvmat.value[l][2] = cv3[l];
-               }
+               //? for(int l=0; l<3; l++) { cvmat.value[l][0] = cv1[l]; cvmat.value[l][1] = cv2[l]; cvmat.value[l][2] = cv3[l]; }
                pstart = (cvmat * pstart);
                pend = (cvmat * pend);
 
@@ -417,7 +442,9 @@ void ParticleTracer::getTriangles( vector<ntlTriangle> *triangles,
                }
        }
 
-       //} // trail
+       } // t
+
+       debMsgStd("ParticleTracer::getTriangles",DM_MSG,"Dumped "<<pcnt<<"/"<<mParts.size()<<" parts, tris:"<<tris<<", showonly:"<<mShowOnly,10);
        return; // DEBUG
 
 #endif // ELBEEM_PLUGIN
index 453c1a7..f08bd58 100644 (file)
@@ -22,10 +22,12 @@ template<class Scalar> class ntlMatrix4x4;
 #define PART_IN     (1<< 8)
 #define PART_OUT    (1<< 9)
 #define PART_INACTIVE (1<<10)
+#define PART_OUTFLUID  (1<<11)
 
 // defines for particle movement
 #define MOVE_FLOATS 1
 #define FLOAT_JITTER 0.03
+//#define FLOAT_JITTER 0.0
 
 extern int ParticleObjectIdCnt;
 
@@ -47,6 +49,8 @@ class ParticleObject
                //! add vector to position
                inline void advance(float vx, float vy, float vz) {
                        mPos[0] += vx; mPos[1] += vy; mPos[2] += vz; }
+               inline void advanceVec(ntlVec3Gfx v) {
+                       mPos[0] += v[0]; mPos[1] += v[1]; mPos[2] += v[2]; }
                //! advance with own velocity
                inline void advanceVel() { mPos += mVel; }
                //! add acceleration to velocity
@@ -83,10 +87,17 @@ class ParticleObject
                        if(set) mStatus &= (~PART_INACTIVE);    
                        else mStatus |= PART_INACTIVE;
                }
+               //! get influid flag
+               inline bool getInFluid() const { return ((mStatus&PART_OUTFLUID)==0); }
+               //! set influid flag
+               inline void setInFluid(bool set) { 
+                       if(set) mStatus &= (~PART_OUTFLUID);    
+                       else mStatus |= PART_OUTFLUID;
+               }
                //! get/set lifetime
-               inline int getLifeTime() const { return mLifeTime; }
+               inline float getLifeTime() const { return mLifeTime; }
                //! set type (lower byte)
-               inline void setLifeTime(int set) { mLifeTime = set; }
+               inline void setLifeTime(float set) { mLifeTime = set; }
                
                inline int getId() const { return mId; }
 
@@ -103,7 +114,7 @@ class ParticleObject
                /*! particle status */
                int mStatus;
                /*! count survived time steps */
-               int mLifeTime;
+               float mLifeTime;
 };
 
 
@@ -161,24 +172,34 @@ class ParticleTracer :
                /*! set/get dump flag */
                inline void setDumpParts(bool set) { mDumpParts = set; }
                inline bool getDumpParts()         { return mDumpParts; }
-               /*! set/get dump flag */
-               inline void setDumpText(bool set) { mDumpText = set; }
-               inline bool getDumpText()         { return mDumpText; }
                /*! set/get dump text file */
-               inline void setDumpTextFile(std::string set) { mDumpTextFile = set; }
-               inline std::string getDumpTextFile()         { return mDumpTextFile; }
+               inline void setDumpTextFile(string set) { mDumpTextFile = set; }
+               inline string getDumpTextFile()         { return mDumpTextFile; }
+               /*! set/get dump text interval */
+               inline void setDumpTextInterval(float set) { mDumpTextInterval = set; }
+               inline float getDumpTextInterval()         { return mDumpTextInterval; }
+               /*! set/get init times */
+               inline void setInitStart(float set) { mInitStart = set; }
+               inline float getInitStart()         { return mInitStart; }
+               inline void setInitEnd(float set)   { mInitEnd = set; }
+               inline float getInitEnd()           { return mInitEnd; }
                
                //! set the particle scaling factor
                inline void setPartScale(float set) { mPartScale = set; }
 
+               //! called after each frame to check if positions should be dumped
+               void checkDumpTextPositions(double simtime);
 
                // NTL geometry implementation
                /*! Get the triangles from this object */
-               virtual void getTriangles( vector<ntlTriangle> *triangles, 
+               virtual void getTriangles(double t, vector<ntlTriangle> *triangles, 
                                vector<ntlVec3Gfx> *vertices, 
                                vector<ntlVec3Gfx> *normals, int objectId );
 
                virtual void notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename,double simtime);
+
+               // notify of next step for trails
+               void checkTrails(double time);
                // free deleted particles
                void cleanup();
 
@@ -211,10 +232,12 @@ class ParticleTracer :
 
                /*! dump particles (or certain types of) to disk? */
                int mDumpParts;
-               /*! dump particles (or certain types of) to disk? */
-               int mDumpText;
                /*! text dump output file */
-               std::string mDumpTextFile;
+               string mDumpTextFile;
+               /*! text dump interval, start at t=0, dumping active if >0 */
+               float mDumpTextInterval;
+               float mDumpTextLastTime;
+               int mDumpTextCount;
                /*! show only a certain type (debugging) */
                int mShowOnly;
                /*! no. of particles to init */
@@ -222,9 +245,17 @@ class ParticleTracer :
 
                //! transform matrix
                ntlMatrix4x4<gfxReal> *mpTrafo;
-
                /*! init sim/pos transformation */
                void initTrafoMatrix();
+
+               //! init time distribution start/end
+               float mInitStart, mInitEnd;
+
+               /*! the particle array (for multiple timesteps) */
+               vector< vector<ParticleObject> > mPrevs;
+               /* prev pos save interval */
+               float mTrailTimeLast, mTrailInterval;
+               int mTrailLength;
 };
 
 #define NTL_PARTICLETRACER_H
index 116d68d..5d8a73e 100644 (file)
@@ -32,11 +32,10 @@ LbmSolverInterface* createSolver();
 SimulationObject::SimulationObject() :
        ntlGeometryShader(),
        mGeoStart(-100.0), mGeoEnd(100.0),
-       mGeoInitId(-1), mpGiTree(NULL), mpGiObjects(NULL),
+       mpGiTree(NULL), mpGiObjects(NULL),
        mpGlob(NULL),
        mPanic( false ),
        mDebugType( 1 /* =FLUIDDISPNothing*/ ),
-       mStepsPerFrame( 10 ),
        mpLbm(NULL), mpParam( NULL ),
        mShowSurface(true), mShowParticles(false),
        mSelectedCid( NULL ),
@@ -48,7 +47,6 @@ SimulationObject::SimulationObject() :
 
        // reset time
        mTime                                           = 0.0;
-       mDisplayTime                    = 0.0;
 }
 
 
@@ -69,12 +67,11 @@ SimulationObject::~SimulationObject()
 /*****************************************************************************/
 /*! init tree for certain geometry init */
 /*****************************************************************************/
-void SimulationObject::initGeoTree(int id) {
+void SimulationObject::initGeoTree() {
        if(mpGlob == NULL) { 
                errFatal("SimulationObject::initGeoTree error","Requires globals!", SIMWORLD_INITERROR); 
                return;
        }
-       mGeoInitId = id;
        ntlScene *scene = mpGlob->getSimScene();
        mpGiObjects = scene->getObjects();
 
@@ -97,6 +94,8 @@ void SimulationObject::freeGeoTree() {
 void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) {
        mpElbeemSettings = new elbeemSimulationSettings;
        *mpElbeemSettings = *settings;
+
+       mGeoInitId = settings->domainId+1;
 }
 
 /******************************************************************************
@@ -104,7 +103,7 @@ void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) {
  *****************************************************************************/
 int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 {
-       if(!SIMWORLD_OK()) return 1;
+       if(! isSimworldOk() ) return 1;
        
        // already inited?
        if(mpLbm) return 0;
@@ -116,16 +115,16 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        }
 
 
+       this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
        //mDimension, mSolverType are deprecated
        string mSolverType(""); 
        mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false ); 
-       //errFatal("SimulationObject::initializeLbmSimulation","Invalid solver type - note that mDimension is deprecated, use the 'solver' keyword instead", SIMWORLD_INITERROR); return 1;
 
        mpLbm = createSolver(); 
   /* check lbm pointer */
        if(mpLbm == NULL) {
                errFatal("SimulationObject::initializeLbmSimulation","Unable to init LBM solver! ", SIMWORLD_INITERROR);
-               return 1;
+               return 2;
        }
        debMsgStd("SimulationObject::initialized",DM_MSG,"IdStr:"<<mpLbm->getIdString() <<" LBM solver! ", 2);
 
@@ -140,10 +139,10 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        mpParts = new ParticleTracer();
        mpParts->parseAttrList( getAttributeList() );
 
-       if(!SIMWORLD_OK()) return 1;
+       if(! isSimworldOk() ) return 3;
        mpParts->setName( getName() + "_part" );
        mpParts->initialize( glob );
-       if(!SIMWORLD_OK()) return 1;
+       if(! isSimworldOk() ) return 4;
        
        // init material settings
        string matMc("default");
@@ -152,6 +151,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        mShowParticles = mpAttrs->readInt("showparticles", mShowParticles, "SimulationObject","mShowParticles", false ); 
 
        checkBoundingBox( mGeoStart, mGeoEnd, "SimulationObject::initializeSimulation" );
+       mpLbm->setLbmInitId( mGeoInitId );
        mpLbm->setGeoStart( mGeoStart );
        mpLbm->setGeoEnd( mGeoEnd );
        mpLbm->setRenderGlobals( mpGlob );
@@ -159,6 +159,8 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        mpLbm->setParticleTracer( mpParts );
        if(mpElbeemSettings) {
                // set further settings from API struct init
+               if(mpElbeemSettings->outputPath) this->mOutFilename = string(mpElbeemSettings->outputPath);
+               mpLbm->initDomainTrafo( mpElbeemSettings->surfaceTrafo );
                mpLbm->setSmoothing(1.0 * mpElbeemSettings->surfaceSmoothing, 1.0 * mpElbeemSettings->surfaceSmoothing);
                mpLbm->setSizeX(mpElbeemSettings->resolutionxyz);
                mpLbm->setSizeY(mpElbeemSettings->resolutionxyz);
@@ -166,11 +168,13 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
                mpLbm->setPreviewSize(mpElbeemSettings->previewresxyz);
                mpLbm->setRefinementDesired(mpElbeemSettings->maxRefine);
                mpLbm->setGenerateParticles(mpElbeemSettings->generateParticles);
+               // set initial particles
+               mpParts->setNumInitialParticles(mpElbeemSettings->numTracerParticles);
 
-               string dinitType = std::string("no");
-               if     (mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_PARTSLIP) dinitType = std::string("part"); 
-               else if(mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_FREESLIP) dinitType = std::string("free"); 
-               else /*if(mpElbeemSettings->obstacleType==FLUIDSIM_OBSTACLE_NOSLIP)*/ dinitType = std::string("no"); 
+               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"); 
                mpLbm->setDomainBound(dinitType);
                mpLbm->setDomainPartSlip(mpElbeemSettings->obstaclePartslip);
                mpLbm->setDumpVelocities(mpElbeemSettings->generateVertexVectors);
@@ -179,7 +183,13 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 
                debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<<mpLbm->getGenerateParticles(),10);
        }
-       mpLbm->initializeSolver();
+
+       if(! mpLbm->initializeSolverMemory()   )         { errMsg("SimulationObject::initialize","initializeSolverMemory failed"); mPanic=true; return 10; }
+       if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverMemory status"); mPanic=true; return 11; } 
+       if(! mpLbm->initializeSolverGrids()    )         { errMsg("SimulationObject::initialize","initializeSolverGrids  failed"); mPanic=true; return 12; }
+       if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverGrids  status"); mPanic=true; return 13; } 
+       if(! mpLbm->initializeSolverPostinit() )         { errMsg("SimulationObject::initialize","initializeSolverPostin failed"); mPanic=true; return 14; }
+       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;
@@ -202,7 +212,6 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        }
        mpLbm->deleteCellIterator( &cid );
 
-#if ELBEEM_PLUGIN!=1
        char charNl = '\n';
        debugOutNnl("SimulationObject::initializeLbmSimulation celltype stats: " <<charNl, 5);
        debugOutNnl("no. of cells = "<<totalCells<<", "<<charNl ,5);
@@ -215,7 +224,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
        }
        // 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];
@@ -226,11 +235,10 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
                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);
+                       debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="<<diffInits<<"!" , 5);
                }
                debugOutNnl(out.str(), 5);
        }
-#endif // ELBEEM_PLUGIN==1
 
        // might be modified by mpLbm
        mpParts->setStart( mGeoStart );
@@ -253,8 +261,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 #ifdef ELBEEM_PLUGIN
        mShowParticles=1;
 #endif // ELBEEM_PLUGIN
-       //if(getenv("ELBEEM_DUMPPARTICLE")) {   // DEBUG ENABLE!!!!!!!!!!
-       if(mpLbm->getGenerateParticles()>0.0) {
+       if((mpLbm->getGenerateParticles()>0.0)||(mpParts->getNumInitialParticles()>0)) {
                mShowParticles=1;
                mpParts->setDumpParts(true);
        }
@@ -293,6 +300,11 @@ void SimulationObject::step( void )
        }
        if(mpLbm->getPanic()) mPanic = true;
 
+       checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0);
+       //if((mpElbeemSettings)&&(mpElbeemSettings->runsimCallback)) {
+               //int ret = (mpElbeemSettings->runsimCallback)(mpElbeemSettings->runsimUserData, FLUIDSIM_CBSTATUS_STEP, 0);
+               //errMsg("runSimulationCallback cbtest1"," "<<this->getName()<<" ret="<<ret);
+       //}
   //debMsgStd("SimulationObject::step",DM_MSG," Sim '"<<mName<<"' stepped to "<<mTime<<" (stept="<<(mpParam->getTimestep())<<", framet="<<getFrameTime()<<") ", 10);
 }
 /*! prepare visualization of simulation for e.g. raytracing */
@@ -399,6 +411,34 @@ void SimulationObject::setMouseClick()
 /*! notify object that dump is in progress (e.g. for field dump) */
 void SimulationObject::notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
        if(!mpLbm) return;
+
        mpLbm->notifySolverOfDump(dumptype, frameNr,frameNrStr,outfilename);
+       checkCallerStatus(FLUIDSIM_CBSTATUS_NEWFRAME, frameNr);
+}
+
+/*! check status (e.g. stop/abort) from calling program, returns !=0 if sth. happened... */
+int SimulationObject::checkCallerStatus(int status, int frame) {
+       //return 0; // DEBUG
+       int ret = 0;
+       if((mpElbeemSettings)&&(mpElbeemSettings->runsimCallback)) {
+               ret = (mpElbeemSettings->runsimCallback)(mpElbeemSettings->runsimUserData, status,frame);
+               if(ret!=FLUIDSIM_CBRET_CONTINUE) {
+                       if(ret==FLUIDSIM_CBRET_STOP) {
+                               debMsgStd("SimulationObject::notifySolverOfDump",DM_NOTIFY,"Got stop signal from caller",1);
+                               setElbeemState( SIMWORLD_STOP );
+                       }
+                       else if(ret==FLUIDSIM_CBRET_ABORT) {
+                               errFatal("SimulationObject::notifySolverOfDump","Got abort signal from caller, aborting...", SIMWORLD_GENERICERROR );
+                               mPanic = 1;
+                       }
+                       else {
+                               errMsg("SimulationObject::notifySolverOfDump","Invalid callback return value: "<<ret<<", ignoring... ");
+                       }
+               }
+       }
+
+       errMsg("SimulationObject::checkCallerStatus","s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
+       if(isSimworldOk()) return 0;
+       return 1;
 }
 
index b8d1e90..200cad8 100644 (file)
@@ -55,7 +55,7 @@ class SimulationObject :
 
 
                /*! init tree for certain geometry init */
-               void initGeoTree(int id);
+               void initGeoTree();
                /*! destroy tree etc. when geometry init done */
                void freeGeoTree();
                /*! get fluid init type at certain position */
@@ -67,10 +67,6 @@ class SimulationObject :
 
                /*! get current (max) simulation time */
                double getCurrentTime( void ) { return mTime; }
-
-               /*! set time to display */
-               void setDisplayTime(double set) { mDisplayTime = set; }
-
                /*! set geometry generation start point */
                virtual void setGeoStart(ntlVec3Gfx set) { mGeoStart = set; }
                /*! set geometry generation end point */
@@ -126,14 +122,11 @@ class SimulationObject :
                /*! current time in the simulation */
                double mTime;
 
-               /*! time to display in the visualizer */
-               double mDisplayTime;
-
                /*! for display - start and end vectors for geometry */
                ntlVec3Gfx mGeoStart, mGeoEnd;
 
                /*! geometry init id */
-               int mGeoInitId;
+               //? int mGeoInitId;
                /*! tree object for geomerty initialization */
                ntlTree *mpGiTree;
                /*! object vector for geo init */
@@ -147,18 +140,9 @@ class SimulationObject :
                /*! debug info to display */
                int mDebugType;
 
-               //! dimension of the simulation - now given by LBM-DIM define globally
-               //! solver type
-
-               /*! when no parametrizer, use this as no. of steps per frame */
-               int mStepsPerFrame;
-
                /*! pointer to the lbm solver */
                LbmSolverInterface *mpLbm;
 
-               /*! marching cubes object */
-               //mCubes *mpMC;
-
                /*! parametrizer for lbm solver */
                Parametrizer *mpParam;
 
@@ -193,8 +177,9 @@ class SimulationObject :
                /*! init parametrizer for anim step length */
                Parametrizer *getParametrizer() { return mpParam; }
 
-               /*! Access marching cubes object */
-               //mCubes *getMCubes( void ) { return mpMC; }
+               /*! check status (e.g. stop/abort) from calling program, returns !=0 if sth. happened... */
+               // parameters same as elbeem runsimCallback
+               int checkCallerStatus(int status, int frame);
 
                /*! get bounding box of fluid for GUI */
                virtual inline ntlVec3Gfx *getBBStart()         { return &mGeoStart; }
diff --git a/intern/elbeem/intern/solver_adap.cpp b/intern/elbeem/intern/solver_adap.cpp
new file mode 100644 (file)
index 0000000..aade9b3
--- /dev/null
@@ -0,0 +1,1293 @@
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * Copyright 2003,2004,2005 Nils Thuerey
+ *
+ * Adaptivity functions
+ *
+ *****************************************************************************/
+
+#include "solver_class.h"
+#include "solver_relax.h"
+#include "particletracer.h"
+
+
+
+/*****************************************************************************/
+//! coarse step functions
+/*****************************************************************************/
+
+
+
+void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
+{
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       lev =0; // get rid of warnings...
+#else
+       FSGR_FORIJK_BOUNDS(lev) {
+               if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
+                       if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
+                               LbmFloat totArea = mFsgrCellArea[0]; // for l=0
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
+                                       if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)&
+                                                       (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
+                                                       ) { 
+                                               totArea += mFsgrCellArea[l];
+                                       }
+                               } // l
+                               QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea;
+                               //continue;
+                       } else
+                       if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & (CFEmpty|CFUnused)) {
+                               QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0;
+                               //continue;
+                       } else {
+                               QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0;
+                       }
+               //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) ); 
+               }
+       } // } TEST DEBUG
+       if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
+#endif //! LBM_NOADCOARSENING==1
+}
+       
+void LbmFsgrSolver::coarseAdvance(int lev)
+{
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       lev =0; // get rid of warnings...
+#else
+       LbmFloat calcCurrentMass = 0.0;
+       LbmFloat calcCurrentVolume = 0.0;
+
+       LbmFloat *ccel = NULL;
+       LbmFloat *tcel = NULL;
+       LbmFloat m[LBM_DFNUM];
+       LbmFloat rho, ux, uy, uz, tmp, usqr, lcsmqo;
+#if OPT3D==1 
+       LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
+#endif // OPT3D==true 
+       m[0] = tmp = usqr = 0.0;
+
+       coarseCalculateFluxareas(lev);
+       // copied from fineAdv.
+       CellFlagType *pFlagSrc = &RFLAG(lev, 1,1,getForZMin1(),SRCS(lev));
+       CellFlagType *pFlagDst = &RFLAG(lev, 1,1,getForZMin1(),TSET(lev));
+       pFlagSrc -= 1;
+       pFlagDst -= 1;
+       ccel = RACPNT(lev, 1,1,getForZMin1() ,SRCS(lev)); // QTEST
+       ccel -= QCELLSTEP;
+       tcel = RACPNT(lev, 1,1,getForZMin1() ,TSET(lev)); // QTEST
+       tcel -= QCELLSTEP;
+       //if(strstr(this->getName().c_str(),"Debug")){ errMsg("DEBUG","DEBUG!!!!!!!!!!!!!!!!!!!!!!!"); }
+
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+#if FSGR_STRICT_DEBUG==1
+               rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
+#endif
+               pFlagSrc++;
+               pFlagDst++;
+               ccel += QCELLSTEP;
+               tcel += QCELLSTEP;
+
+               // from coarse cells without unused nbs are not necessary...! -> remove
+               if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
+                       bool invNb = false;
+                       FORDF1 { if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; } }   
+                       if(!invNb) {
+                               // WARNING - this modifies source flag array...
+                               *pFlagSrc = CFFluid|CFGrNorm;
+#if ELBEEM_PLUGIN!=1
+                               errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
+#endif // ELBEEM_PLUGIN!=1
+                               // move to perform coarsening?
+                       }
+               } // */
+
+#if FSGR_STRICT_DEBUG==1
+               *pFlagDst = *pFlagSrc; // always set other set...
+#else
+               *pFlagDst = (*pFlagSrc & (~CFGrCoarseInited)); // always set other set... , remove coarse inited flag
+#endif
+
+               // old INTCFCOARSETEST==1
+               if((*pFlagSrc) & CFGrFromCoarse) {  // interpolateFineFromCoarse test!
+                       if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
+                               FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
+                       } else {
+                               interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
+                               this->mNumUsedCells++;
+                       }
+                       continue; // interpolateFineFromCoarse test!
+               } // interpolateFineFromCoarse test! old INTCFCOARSETEST==1
+
+               if( ((*pFlagSrc) & (CFFluid)) ) { 
+                       ccel = RACPNT(lev, i,j,k ,SRCS(lev)); 
+                       tcel = RACPNT(lev, i,j,k ,TSET(lev));
+
+                       if( ((*pFlagSrc) & (CFGrFromFine)) ) { 
+                               FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
+                               continue; // comes from fine grid
+                       }
+                       // also ignore CFGrFromCoarse
+                       else if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
+                               FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
+                               continue; 
+                       }
+
+                       OPTIMIZED_STREAMCOLLIDE;
+                       *pFlagDst |= CFNoBndFluid; // test?
+                       calcCurrentVolume += RAC(ccel,dFlux); 
+                       calcCurrentMass   += RAC(ccel,dFlux)*rho;
+                       //ebugMarkCell(lev+1, 2*i+1,2*j+1,2*k  );
+#if FSGR_STRICT_DEBUG==1
+                       if(rho<-1.0){ debugMarkCell(lev, i,j,k ); 
+                               errMsg("INVRHOCELL_CHECK"," l"<<lev<<" "<< PRINT_IJK<<" rho:"<<rho ); 
+                               CAUSE_PANIC;
+                       }
+#endif // FSGR_STRICT_DEBUG==1
+                       this->mNumUsedCells++;
+
+               }
+       } 
+       pFlagSrc+=2; // after x
+       pFlagDst+=2; // after x
+       ccel += (QCELLSTEP*2);
+       tcel += (QCELLSTEP*2);
+       } 
+       pFlagSrc+= mLevel[lev].lSizex*2; // after y
+       pFlagDst+= mLevel[lev].lSizex*2; // after y
+       ccel += (QCELLSTEP*mLevel[lev].lSizex*2);
+       tcel += (QCELLSTEP*mLevel[lev].lSizex*2);
+       } // all cell loop k,j,i
+       
+
+       //errMsg("coarseAdvance","level "<<lev<<" stepped from "<<mLevel[lev].setCurr<<" to "<<mLevel[lev].setOther);
+       if(!this->mSilent){ errMsg("coarseAdvance","level "<<lev<<" stepped from "<<SRCS(lev)<<" to "<<TSET(lev)); }
+       // */
+
+       // update other set
+  mLevel[lev].setOther   = mLevel[lev].setCurr;
+  mLevel[lev].setCurr   ^= 1;
+  mLevel[lev].lsteps++;
+  mLevel[lev].lmass   = calcCurrentMass   * mLevel[lev].lcellfactor;
+  mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
+#if ELBEEM_PLUGIN!=1
+  errMsg("DFINI", " m l"<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor );
+  errMsg("DFINI", " v l"<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor );
+#endif // ELBEEM_PLUGIN
+#endif //! LBM_NOADCOARSENING==1
+}
+
+/*****************************************************************************/
+//! multi level functions
+/*****************************************************************************/
+
+
+// get dfs from level (lev+1) to (lev) coarse border nodes
+void LbmFsgrSolver::coarseRestrictFromFine(int lev)
+{
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       lev =0; // get rid of warnings...
+#else
+       if((lev<0) || ((lev+1)>mMaxRefine)) return;
+#if FSGR_STRICT_DEBUG==1
+       // reset all unused cell values to invalid
+       int unuCnt = 0;
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+       for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+       for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+               CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,mLevel[lev].setCurr);
+               if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
+                       FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
+                       unuCnt++;
+                       // set here
+               } else if( ((*pFlagSrc) & (CFFluid|CFGrNorm)) == (CFFluid|CFGrNorm) ) { 
+                       // simulated...
+               } else {
+                       // reset in interpolation
+                       //errMsg("coarseRestrictFromFine"," reset l"<<lev<<" "<<PRINT_IJK);
+               }
+               if( ((*pFlagSrc) & (CFEmpty|CFUnused)) ) {  // test, also reset?
+                       FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
+               } // test
+       } } }
+       errMsg("coarseRestrictFromFine"," reset l"<<lev<<" fluid|coarseBorder cells: "<<unuCnt);
+#endif // FSGR_STRICT_DEBUG==1
+       const int srcSet = mLevel[lev+1].setCurr;
+       const int dstSet = mLevel[lev].setCurr;
+
+       //restrict
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+       for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+       for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+               CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,dstSet);
+               if((*pFlagSrc) & (CFFluid)) { 
+                       if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
+                               // do resctriction
+                               mNumInterdCells++;
+                               coarseRestrictCell(lev, i,j,k,srcSet,dstSet);
+
+                               this->mNumUsedCells++;
+                       } // from fine & fluid
+                       else {
+                               if(RFLAG(lev+1, 2*i,2*j,2*k,srcSet) & CFGrFromCoarse) {
+                                       RFLAG(lev, i,j,k,dstSet) |= CFGrToFine;
+                               } else {
+                                       RFLAG(lev, i,j,k,dstSet) &= (~CFGrToFine);
+                               }
+                       }
+               } // & fluid
+       }}}
+       if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
+#endif //! LBM_NOADCOARSENING==1
+}
+
+bool LbmFsgrSolver::adaptGrid(int lev) {
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       lev =0; // get rid of warnings...
+       return false;
+#else
+       if((lev<0) || ((lev+1)>mMaxRefine)) return false;
+       bool change = false;
+       { // refinement, PASS 1-3
+
+       //bool nbsok;
+       // FIXME remove TIMEINTORDER ?
+       LbmFloat interTime = 0.0;
+       // update curr from other, as streaming afterwards works on curr
+       // thus read only from srcSet, modify other
+       const int srcSet = mLevel[lev].setOther;
+       const int dstSet = mLevel[lev].setCurr;
+       const int srcFineSet = mLevel[lev+1].setCurr;
+       const bool debugRefinement = false;
+
+       // use //template functions for 2D/3D
+                       /*if(strstr(this->getName().c_str(),"Debug"))
+                       if(lev+1==mMaxRefine) { // mixborder
+                               for(int l=0;((l<this->cDirNum) && (!removeFromFine)); l++) {  // FARBORD
+                                       int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
+                                       if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&CFBnd) { // NEWREFT
+                                               removeFromFine=true;
+                                       }
+                               }
+                       } // FARBORD */
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+
+               if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
+                       bool removeFromFine = false;
+                       const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
+                       CellFlagType reqType = CFGrNorm;
+                       if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
+                       
+                       if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & reqType) &&
+                           (!(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & (notAllowed)) )  ){ // ok
+                       } else {
+                               removeFromFine=true;
+                       }
+
+                       if(removeFromFine) {
+                               // dont turn CFGrFromFine above interface cells into CFGrNorm
+                               //errMsg("performRefinement","Removing CFGrFromFine on lev"<<lev<<" " <<PRINT_IJK<<" srcflag:"<<convertCellFlagType2String(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet)) <<" set:"<<dstSet );
+                               RFLAG(lev, i,j,k, dstSet) = CFEmpty;
+#if FSGR_STRICT_DEBUG==1
+                               // for interpolation later on during fine grid fixing
+                               // these cells are still correctly inited
+                               RFLAG(lev, i,j,k, dstSet) |= CFGrCoarseInited;  // remove later on? FIXME?
+#endif // FSGR_STRICT_DEBUG==1
+                               //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFEmpty; // FLAGTEST
+                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,i,j,k); 
+                               change=true;
+                               mNumFsgrChanges++;
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+                                       //errMsg("performRefinement","On lev:"<<lev<<" check: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet<<" = "<<convertCellFlagType2String(RFLAG(lev, ni,nj,nk, srcSet)) );
+                                       if( (  RFLAG(lev, ni,nj,nk, srcSet)&CFFluid      ) &&
+                                                       (!(RFLAG(lev, ni,nj,nk, srcSet)&CFGrFromFine)) ) { // dont change status of nb. from fine cells
+                                               // tag as inited for debugging, cell contains fluid DFs anyway
+                                               RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromFine|CFGrCoarseInited;
+                                               //errMsg("performRefinement","On lev:"<<lev<<" set to from fine: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet);
+                                               //if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
+                                       }
+                               } // l 
+
+                               // FIXME fix fine level?
+                       }
+
+                       // recheck from fine flag
+               }
+       }}} // TEST
+       // PASS 1 */
+
+
+               /*if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
+                       bool invNb = false;
+                       FORDF1 { 
+                               if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; }
+                       }   
+                       if(!invNb) {
+                               *pFlagSrc = CFFluid|CFGrNorm;
+                               errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
+                       }
+               } // */
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
+
+               // test from coarseAdvance
+               // from coarse cells without unused nbs are not necessary...! -> remove
+
+               if(RFLAG(lev, i,j,k, srcSet) & CFGrFromCoarse) {
+
+                       // from coarse cells without unused nbs are not necessary...! -> remove
+                       bool invNb = false;
+                       bool fluidNb = false;
+                       for(int l=1; l<this->cDirNum; l++) { 
+                               if(RFLAG_NB(lev, i, j, k, srcSet, l) & CFUnused) { invNb = true; }
+                               if(RFLAG_NB(lev, i, j, k, srcSet, l) & (CFGrNorm)) { fluidNb = true; }
+                       }   
+                       if(!invNb) {
+                               // no unused cells around -> calculate normally from now on
+                               RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
+                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
+                               change=true;
+                               mNumFsgrChanges++;
+                       } // from advance 
+                       if(!fluidNb) {
+                               // no fluid cells near -> no transfer necessary
+                               RFLAG(lev, i,j,k, dstSet) = CFUnused;
+                               //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFUnused; // FLAGTEST
+                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
+                               change=true;
+                               mNumFsgrChanges++;
+                       } // from advance 
+
+
+                       // dont allow double transfer
+                       // this might require fixing the neighborhood
+                       if(RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse)) { 
+                               // dont turn CFGrFromFine above interface cells into CFGrNorm
+                               //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<lev<<" " <<PRINT_IJK<<" due to finer from coarse cell " );
+                               RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
+                               if(lev>0) RFLAG(lev-1, i/2,j/2,k/2, mLevel[lev-1].setCurr) &= (~CFGrToFine); // TODO add more of these?
+                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
+                               change=true;
+                               mNumFsgrChanges++;
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+                                       if(RFLAG(lev, ni,nj,nk, srcSet)&(CFGrNorm)) { //ok
+                                               for(int m=1; m<this->cDirNum; m++) { 
+                                                       int mi=  ni +this->dfVecX[m], mj=  nj +this->dfVecY[m], mk=  nk +this->dfVecZ[m];
+                                                       if(RFLAG(lev,  mi, mj, mk, srcSet)&CFUnused) {
+                                                               // norm cells in neighborhood with unused nbs have to be new border...
+                                                               RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromCoarse;
+                                                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
+                                                       }
+                                               }
+                                               // these alreay have valid values...
+                                       }
+                                       else if(RFLAG(lev, ni,nj,nk, srcSet)&(CFUnused)) { //ok
+                                               // this should work because we have a valid neighborhood here for now
+                                               interpolateCellFromCoarse(lev,  ni, nj, nk, dstSet, interTime, CFFluid|CFGrFromCoarse, false);
+                                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
+                                               mNumFsgrChanges++;
+                                       }
+                               } // l 
+                       } // double transer
+
+               } // from coarse
+
+       } } }
+       // PASS 2 */
+
+
+       // fix dstSet from fine cells here
+       // warning - checks CFGrFromFine on dstset changed before!
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
+
+               //if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
+               if(RFLAG(lev, i,j,k, dstSet) & CFGrFromFine) {
+                       // modify finer level border
+                       if((RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse))) { 
+                               //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" from l"<<lev<<" " <<PRINT_IJK );
+                               CellFlagType setf = CFFluid;
+                               if(lev+1 < mMaxRefine) setf = CFFluid|CFGrNorm;
+                               RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)=setf;
+                               change=true;
+                               mNumFsgrChanges++;
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l];
+                                       if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFGrFromCoarse)) {
+                                               //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
+                                               RFLAG(lev+1,  bi, bj, bk, srcFineSet) = setf;
+                                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
+                                       }
+                                       else if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFUnused      )) { 
+                                               //errMsg("performRefinement","Removing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
+                                               interpolateCellFromCoarse(lev+1,  bi, bj, bk, srcFineSet, interTime, setf, false);
+                                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
+                                               mNumFsgrChanges++;
+                                       }
+                               }
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l];
+                                       if(   (RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFFluid       ) &&
+                                                       (!(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFGrFromCoarse)) ) {
+                                               // all unused nbs now of coarse have to be from coarse
+                                               for(int m=1; m<this->cDirNum; m++) { 
+                                                       int mi=  bi +this->dfVecX[m], mj=  bj +this->dfVecY[m], mk=  bk +this->dfVecZ[m];
+                                                       if(RFLAG(lev+1,  mi, mj, mk, srcFineSet)&CFUnused) {
+                                                               //errMsg("performRefinement","Changing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(mi,mj,mk) );
+                                                               interpolateCellFromCoarse(lev+1,  mi, mj, mk, srcFineSet, interTime, CFFluid|CFGrFromCoarse, false);
+                                                               if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,mi,mj,mk); 
+                                                               mNumFsgrChanges++;
+                                                       }
+                                               }
+                                               // nbs prepared...
+                                       }
+                               }
+                       }
+                       
+               } // convert regions of from fine
+       }}} // TEST
+       // PASS 3 */
+
+       if(!this->mSilent){ errMsg("performRefinement"," for l"<<lev<<" done ("<<change<<") " ); }
+       } // PASS 1-3
+       // refinement done
+
+       //LbmFsgrSolver::performCoarsening(int lev) {
+       { // PASS 4,5
+       bool nbsok;
+       // WARNING
+       // now work on modified curr set
+       const int srcSet = mLevel[lev].setCurr;
+       const int dstlev = lev+1;
+       const int dstFineSet = mLevel[dstlev].setCurr;
+       const bool debugCoarsening = false;
+
+       // PASS 5 test DEBUG
+       /*if(this->mInitDone) {
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+                       if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) {
+                               // check empty -> from fine conversion
+                               bool changeToFromFine = false;
+                               const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
+                               CellFlagType reqType = CFGrNorm;
+                               if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
+                               if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
+                                   (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
+                                       changeToFromFine=true; }
+                               if(changeToFromFine) {
+                                       change = true;
+                                       mNumFsgrChanges++;
+                                       RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine;
+                                       if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
+                                       // same as restr from fine func! not necessary ?!
+                                       // coarseRestrictFromFine part 
+                                       coarseRestrictCell(lev, i,j,k,srcSet, dstFineSet);
+                               }
+                       } // only check empty cells
+       }}} // TEST!
+       } // PASS 5 */
+
+       // use //template functions for 2D/3D
+                                       /*if(strstr(this->getName().c_str(),"Debug"))
+                                       if((nbsok)&&(lev+1==mMaxRefine)) { // mixborder
+                                               for(int l=0;((l<this->cDirNum) && (nbsok)); l++) {  // FARBORD
+                                                       int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
+                                                       if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
+                                                               nbsok=false;
+                                                       }
+                                               }
+                                       } // FARBORD */
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+
+                       // from coarse cells without unused nbs are not necessary...! -> remove
+                       // perform check from coarseAdvance here?
+                       if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
+                               // remove from fine cells now that are completely in fluid
+                               // FIXME? check that new from fine in performRefinement never get deleted here afterwards?
+                               // or more general, one cell never changed more than once?
+                               const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
+                               //const CellFlagType notNbAllowed = (CFInter|CFBnd|CFGrFromFine); unused
+                               CellFlagType reqType = CFGrNorm;
+                               if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
+
+                               nbsok = true;
+                               for(int l=0; l<this->cDirNum && nbsok; l++) { 
+                                       int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
+                                       if(   (RFLAG(lev+1, ni,nj,nk, dstFineSet) & reqType) &&
+                                                       (!(RFLAG(lev+1, ni,nj,nk, dstFineSet) & (notAllowed)) )  ){
+                                               // ok
+                                       } else {
+                                               nbsok=false;
+                                       }
+                                       // FARBORD
+                               }
+                               // dont turn CFGrFromFine above interface cells into CFGrNorm
+                               // now check nbs on same level
+                               for(int l=1; l<this->cDirNum && nbsok; l++) { 
+                                       int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+                                       if(RFLAG(lev, ni,nj,nk, srcSet)&(CFFluid)) { //ok
+                                       } else {
+                                               nbsok = false;
+                                       }
+                               } // l
+
+                               if(nbsok) {
+                                       // conversion to coarse fluid cell
+                                       change = true;
+                                       mNumFsgrChanges++;
+                                       RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrNorm;
+                                       // dfs are already ok...
+                                       //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine changed to CFGrNorm at lev"<<lev<<" " <<PRINT_IJK );
+                                       if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
+
+                                       // only check complete cubes
+                                       for(int dx=-1;dx<=1;dx+=2) {
+                                       for(int dy=-1;dy<=1;dy+=2) {
+                                       for(int dz=-1*(LBMDIM&1);dz<=1*(LBMDIM&1);dz+=2) { // 2d/3d
+                                               // check for norm and from coarse, as the coarse level might just have been refined...
+                                               if( 
+                                                               // we now the flag of the current cell! ( RFLAG(lev, i   , j   , k   ,  srcSet)&(CFGrNorm)) &&
+                                                               ( RFLAG(lev, i+dx, j   , k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+                                                               ( RFLAG(lev, i   , j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+                                                               ( RFLAG(lev, i   , j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+
+                                                               ( RFLAG(lev, i+dx, j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+                                                               ( RFLAG(lev, i+dx, j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+                                                               ( RFLAG(lev, i   , j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
+                                                               ( RFLAG(lev, i+dx, j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) 
+                                                       ) {
+                                                       // middle source node on higher level
+                                                       int dstx = (2*i)+dx;
+                                                       int dsty = (2*j)+dy;
+                                                       int dstz = (2*k)+dz;
+
+                                                       mNumFsgrChanges++;
+                                                       RFLAG(dstlev, dstx,dsty,dstz, dstFineSet) = CFUnused;
+                                                       RFLAG(dstlev, dstx,dsty,dstz, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
+                                                       //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init center unused set l"<<dstlev<<" at "<<PRINT_VEC(dstx,dsty,dstz) );
+
+                                                       for(int l=1; l<this->cDirNum; l++) { 
+                                                               int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l];
+                                                               if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFFluid)) { 
+                                                                       RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
+                                                               }
+                                                               if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFInter)) { 
+                                                                       //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init CHECK Warning - deleting interface cell...");
+                                                                       this->mFixMass += QCELL( dstlev, dstni,dstnj,dstnk, dstFineSet, dMass);
+                                                                       RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
+                                                               }
+                                                       } // l
+
+                                                       // again check nb flags of all surrounding cells to see if any from coarse
+                                                       // can be convted to unused
+                                                       for(int l=1; l<this->cDirNum; l++) { 
+                                                               int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l];
+                                                               // have to be at least from coarse here...
+                                                               //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" "<< convertCellFlagType2String(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)) );
+                                                               if(!(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFUnused) )) { 
+                                                                       bool delok = true;
+                                                                       // careful long range here... check domain bounds?
+                                                                       for(int m=1; m<this->cDirNum; m++) {                                                                            
+                                                                               int chkni=dstni+this->dfVecX[m], chknj=dstnj+this->dfVecY[m], chknk=dstnk+this->dfVecZ[m];
+                                                                               if(RFLAG(dstlev, chkni,chknj,chknk, dstFineSet)&(CFUnused|CFGrFromCoarse)) { 
+                                                                                       // this nb cell is ok for deletion
+                                                                               } else { 
+                                                                                       delok=false; // keep it!
+                                                                               }
+                                                                               //errMsg("performCoarsening"," CHECK "<<PRINT_VEC(dstni,dstnj,dstnk)<<" to "<<PRINT_VEC( chkni,chknj,chknk )<<" f:"<< convertCellFlagType2String( RFLAG(dstlev, chkni,chknj,chknk, dstFineSet))<<" nbsok"<<delok );
+                                                                       }
+                                                                       //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" ok"<<delok );
+                                                                       if(delok) {
+                                                                               mNumFsgrChanges++;
+                                                                               RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFUnused;
+                                                                               RFLAG(dstlev, dstni,dstnj,dstnk, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
+                                                                               if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(dstlev,dstni,dstnj,dstnk); 
+                                                                       }
+                                                               }
+                                                       } // l
+                                                       // treat subcube
+                                                       //ebugMarkCell(lev,i+dx,j+dy,k+dz); 
+                                                       //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init, dir:"<<PRINT_VEC(dx,dy,dz) );
+                                               }
+                                       } } }
+
+                               }   // ?
+                       } // convert regions of from fine
+       }}} // TEST!
+       // PASS 4 */
+
+               // reinit cell area value
+               /*if( RFLAG(lev, i,j,k,srcSet) & CFFluid) {
+                       if( RFLAG(lev+1, i*2,j*2,k*2,dstFineSet) & CFGrFromCoarse) {
+                               LbmFloat totArea = mFsgrCellArea[0]; // for l=0
+                               for(int l=1; l<this->cDirNum; l++) { 
+                                       int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
+                                       if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&
+                                                       (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
+                                                       //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
+                                                       ) { 
+                                               //LbmFloat area = 0.25; if(this->dfVecX[l]!=0) area *= 0.5; if(this->dfVecY[l]!=0) area *= 0.5; if(this->dfVecZ[l]!=0) area *= 0.5;
+                                               totArea += mFsgrCellArea[l];
+                                       }
+                               } // l
+                               QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
+                               QCELL(lev, i,j,k,srcSet, dFlux) = totArea;
+                       } else {
+                               QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
+                               QCELL(lev, i,j,k,srcSet, dFlux) = 1.0;
+                       }
+                       //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,srcSet, dFlux) );
+               }
+       // */
+
+
+       // PASS 5 org
+       /*if(strstr(this->getName().c_str(),"Debug"))
+       if((changeToFromFine)&&(lev+1==mMaxRefine)) { // mixborder
+               for(int l=0;((l<this->cDirNum) && (changeToFromFine)); l++) {  // FARBORD
+                       int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
+                       if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
+                               changeToFromFine=false; }
+               } 
+       }// FARBORD */
+       //if(!this->mInitDone) {
+       for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
+  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
+  for(int i=1;i<mLevel[lev].lSizex-1;++i) {
+
+
+                       if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) {
+                               // check empty -> from fine conversion
+                               bool changeToFromFine = false;
+                               const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
+                               CellFlagType reqType = CFGrNorm;
+                               if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
+
+                               if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
+                                   (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
+                                       // DEBUG 
+                                       changeToFromFine=true;
+                               }
+
+                               // FARBORD
+
+                               if(changeToFromFine) {
+                                       change = true;
+                                       mNumFsgrChanges++;
+                                       RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine;
+                                       if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
+                                       // same as restr from fine func! not necessary ?!
+                                       // coarseRestrictFromFine part 
+                               }
+                       } // only check empty cells
+
+       }}} // TEST!
+       //} // init done
+       // PASS 5 */
+       } // coarsening, PASS 4,5
+
+       if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); }
+       return change;
+#endif //! LBM_NOADCOARSENING==1
+}
+
+/*****************************************************************************/
+//! cell restriction and prolongation
+/*****************************************************************************/
+
+void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet)
+{
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       i=j=k=srcSet=dstSet=lev =0; // get rid of warnings...
+#else
+       LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
+       LbmFloat *tcel = RACPNT(lev  , i,j,k      ,dstSet);
+
+       LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;                       
+       //LbmFloat *ccel = NULL;
+       //LbmFloat *tcel = NULL;
+#if OPT3D==1 
+       LbmFloat m[LBM_DFNUM];
+       // for macro add
+       LbmFloat usqr;
+       //LbmFloat *addfcel, *dstcell;
+       LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
+       LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
+#else // OPT3D==true 
+       LbmFloat df[LBM_DFNUM];
+       LbmFloat omegaDst, omegaSrc;
+       LbmFloat feq[LBM_DFNUM];
+       LbmFloat dfScale = mDfScaleUp;
+#endif // OPT3D==true 
+
+#                              if OPT3D==0
+       // add up weighted dfs
+       FORDF0{ df[l] = 0.0;}
+       for(int n=0;(n<this->cDirNum); n++) { 
+               int ni=2*i+1*this->dfVecX[n], nj=2*j+1*this->dfVecY[n], nk=2*k+1*this->dfVecZ[n];
+               ccel = RACPNT(lev+1, ni,nj,nk,srcSet);// CFINTTEST
+               const LbmFloat weight = mGaussw[n];
+               FORDF0{
+                       LbmFloat cdf = weight * RAC(ccel,l);
+#                                              if FSGR_STRICT_DEBUG==1
+                       if( cdf<-1.0 ){ errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]); }
+#                                              endif
+                       //errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]<<" = "<<cdf<<" , w"<<weight); 
+                       df[l] += cdf;
+               }
+       }
+
+       // calc rho etc. from weighted dfs
+       rho = ux  = uy  = uz  = 0.0;
+       FORDF0{
+               LbmFloat cdf = df[l];
+               rho += cdf; 
+               ux  += (this->dfDvecX[l]*cdf); 
+               uy  += (this->dfDvecY[l]*cdf);  
+               uz  += (this->dfDvecZ[l]*cdf);  
+       }
+
+       FORDF0{ feq[l] = this->getCollideEq(l, rho,ux,uy,uz); }
+       if(mLevel[lev  ].lcsmago>0.0) {
+               const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feq);
+               omegaDst  = this->getLesOmega(mLevel[lev  ].omega,mLevel[lev  ].lcsmago,Qo);
+               omegaSrc = this->getLesOmega(mLevel[lev+1].omega,mLevel[lev+1].lcsmago,Qo);
+       } else {
+               omegaDst = mLevel[lev+0].omega; /* NEWSMAGOT*/ 
+               omegaSrc = mLevel[lev+1].omega;
+       }
+       dfScale   = (mLevel[lev  ].timestep/mLevel[lev+1].timestep)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0); // yu
+       FORDF0{
+               RAC(tcel, l) = feq[l]+ (df[l]-feq[l])*dfScale;
+       } 
+#                              else // OPT3D
+       // similar to OPTIMIZED_STREAMCOLLIDE_UNUSED
+                                                               
+       //rho = ux = uy = uz = 0.0;
+       MSRC_C  = CCELG_C(0) ;
+       MSRC_N  = CCELG_N(0) ;
+       MSRC_S  = CCELG_S(0) ;
+       MSRC_E  = CCELG_E(0) ;
+       MSRC_W  = CCELG_W(0) ;
+       MSRC_T  = CCELG_T(0) ;
+       MSRC_B  = CCELG_B(0) ;
+       MSRC_NE = CCELG_NE(0);
+       MSRC_NW = CCELG_NW(0);
+       MSRC_SE = CCELG_SE(0);
+       MSRC_SW = CCELG_SW(0);
+       MSRC_NT = CCELG_NT(0);
+       MSRC_NB = CCELG_NB(0);
+       MSRC_ST = CCELG_ST(0);
+       MSRC_SB = CCELG_SB(0);
+       MSRC_ET = CCELG_ET(0);
+       MSRC_EB = CCELG_EB(0);
+       MSRC_WT = CCELG_WT(0);
+       MSRC_WB = CCELG_WB(0);
+       for(int n=1;(n<this->cDirNum); n++) { 
+               ccel = RACPNT(lev+1,  2*i+1*this->dfVecX[n], 2*j+1*this->dfVecY[n], 2*k+1*this->dfVecZ[n]  ,srcSet);
+               MSRC_C  += CCELG_C(n) ;
+               MSRC_N  += CCELG_N(n) ;
+               MSRC_S  += CCELG_S(n) ;
+               MSRC_E  += CCELG_E(n) ;
+               MSRC_W  += CCELG_W(n) ;
+               MSRC_T  += CCELG_T(n) ;
+               MSRC_B  += CCELG_B(n) ;
+               MSRC_NE += CCELG_NE(n);
+               MSRC_NW += CCELG_NW(n);
+               MSRC_SE += CCELG_SE(n);
+               MSRC_SW += CCELG_SW(n);
+               MSRC_NT += CCELG_NT(n);
+               MSRC_NB += CCELG_NB(n);
+               MSRC_ST += CCELG_ST(n);
+               MSRC_SB += CCELG_SB(n);
+               MSRC_ET += CCELG_ET(n);
+               MSRC_EB += CCELG_EB(n);
+               MSRC_WT += CCELG_WT(n);
+               MSRC_WB += CCELG_WB(n);
+       }
+       rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  
+               + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT 
+               + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; 
+       ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW 
+               + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB;  
+       uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW 
+               + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB;  
+       uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB 
+               + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB;  
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
+       \
+       lcsmeq[dC] = EQC ; \
+       COLL_CALCULATE_DFEQ(lcsmeq); \
+       COLL_CALCULATE_NONEQTENSOR(lev+0, MSRC_ )\
+       COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \
+       COLL_CALCULATE_CSMOMEGAVAL(lev+1, lcsmSrcOmega); \
+       \
+       lcsmdfscale   = (mLevel[lev+0].timestep/mLevel[lev+1].timestep)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0);  \
+       RAC(tcel, dC ) = (lcsmeq[dC ] + (MSRC_C -lcsmeq[dC ] )*lcsmdfscale);
+       RAC(tcel, dN ) = (lcsmeq[dN ] + (MSRC_N -lcsmeq[dN ] )*lcsmdfscale);
+       RAC(tcel, dS ) = (lcsmeq[dS ] + (MSRC_S -lcsmeq[dS ] )*lcsmdfscale);
+       RAC(tcel, dE ) = (lcsmeq[dE ] + (MSRC_E -lcsmeq[dE ] )*lcsmdfscale);
+       RAC(tcel, dW ) = (lcsmeq[dW ] + (MSRC_W -lcsmeq[dW ] )*lcsmdfscale);
+       RAC(tcel, dT ) = (lcsmeq[dT ] + (MSRC_T -lcsmeq[dT ] )*lcsmdfscale);
+       RAC(tcel, dB ) = (lcsmeq[dB ] + (MSRC_B -lcsmeq[dB ] )*lcsmdfscale);
+       RAC(tcel, dNE) = (lcsmeq[dNE] + (MSRC_NE-lcsmeq[dNE] )*lcsmdfscale);
+       RAC(tcel, dNW) = (lcsmeq[dNW] + (MSRC_NW-lcsmeq[dNW] )*lcsmdfscale);
+       RAC(tcel, dSE) = (lcsmeq[dSE] + (MSRC_SE-lcsmeq[dSE] )*lcsmdfscale);
+       RAC(tcel, dSW) = (lcsmeq[dSW] + (MSRC_SW-lcsmeq[dSW] )*lcsmdfscale);
+       RAC(tcel, dNT) = (lcsmeq[dNT] + (MSRC_NT-lcsmeq[dNT] )*lcsmdfscale);
+       RAC(tcel, dNB) = (lcsmeq[dNB] + (MSRC_NB-lcsmeq[dNB] )*lcsmdfscale);
+       RAC(tcel, dST) = (lcsmeq[dST] + (MSRC_ST-lcsmeq[dST] )*lcsmdfscale);
+       RAC(tcel, dSB) = (lcsmeq[dSB] + (MSRC_SB-lcsmeq[dSB] )*lcsmdfscale);
+       RAC(tcel, dET) = (lcsmeq[dET] + (MSRC_ET-lcsmeq[dET] )*lcsmdfscale);
+       RAC(tcel, dEB) = (lcsmeq[dEB] + (MSRC_EB-lcsmeq[dEB] )*lcsmdfscale);
+       RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
+       RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
+#                              endif // OPT3D==0
+#endif //! LBM_NOADCOARSENING==1
+}
+
+void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
+#if LBM_NOADCOARSENING==1
+       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
+       i=j=k=dstSet=lev =0; // get rid of warnings...
+       t=0.0; flagSet=0; markNbs=false;
+#else
+       LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
+       LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+
+#if OPT3D==1 
+       // for macro add
+       LbmFloat addDfFacT, addVal, usqr;
+       LbmFloat *addfcel, *dstcell;
+       LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
+       LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
+#endif // OPT3D==true 
+
+       // SET required nbs to from coarse (this might overwrite flag several times)
+       // this is not necessary for interpolateFineFromCoarse
+       if(markNbs) {
+       FORDF1{ 
+               int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
+               if(RFLAG(lev,ni,nj,nk,dstSet)&CFUnused) {
+                       // parents have to be inited!
+                       interpolateCellFromCoarse(lev, ni, nj, nk, dstSet, t, CFFluid|CFGrFromCoarse, false);
+               }
+       } }
+
+       // change flag of cell to be interpolated
+       RFLAG(lev,i,j,k, dstSet) = flagSet;
+       mNumInterdCells++;
+
+       // interpolation lines...
+       int betx = i&1;
+       int bety = j&1;
+       int betz = k&1;
+       
+       if((!betx) && (!bety) && (!betz)) {
+               ADD_INT_DFS(lev-1, i/2  ,j/2  ,k/2  , 0.0, 1.0);
+       }
+       else if(( betx) && (!bety) && (!betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D1);
+       }
+       else if((!betx) && ( bety) && (!betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D1);
+       }
+       else if((!betx) && (!bety) && ( betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D1);
+       }
+       else if(( betx) && ( bety) && (!betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D2);
+       }
+       else if((!betx) && ( bety) && ( betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D2);
+       }
+       else if(( betx) && (!bety) && ( betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D2);
+       }
+       else if(( betx) && ( bety) && ( betz)) {
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D3);
+               ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)+1, t, WO1D3);
+       }
+       else {
+               CAUSE_PANIC;
+               errFatal("interpolateCellFromCoarse","Invalid!?", SIMWORLD_GENERICERROR);       
+       }
+
+       IDF_WRITEBACK;
+       return;
+#endif //! LBM_NOADCOARSENING==1
+}
+
+
+
+/*****************************************************************************/
+/*! change the  size of the LBM time step */
+/*****************************************************************************/
+void LbmFsgrSolver::adaptTimestep() {
+       LbmFloat massTOld=0.0, massTNew=0.0;
+       LbmFloat volTOld=0.0, volTNew=0.0;
+
+       bool rescale = false;  // do any rescale at all?
+       LbmFloat scaleFac = -1.0; // timestep scaling
+       if(this->mPanic) return;
+
+       LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS];
+       LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS];
+       for(int lev=mMaxRefine; lev>=0 ; lev--) {
+               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);
+
+       //newdt = this->mpParam->getTimestep() * (allowMax/nextmax);
+       LbmFloat newdt = this->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;
+       } else {
+               if(nextmax<allowMax*reduceFac) {
+                       newdt = this->mpParam->getTimestep() / reduceFac;
+               }
+       } // newtr
+       //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< this->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(nextmax>allowMax/reduceFac){ minCutoff=true; } // only if really large vels...
+       }
+
+       LbmFloat dtdiff = fabs(newdt - this->mpParam->getTimestep());
+       if(!this->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); }
+
+       // in range, and more than X% change?
+       //if( newdt <  this->mpParam->getTimestep() ) // DEBUG
+       LbmFloat rhoAvg = mCurrentMass/mCurrentVolume;
+       if( (newdt<=this->mpParam->getMaxTimestep()) && (newdt>=this->mpParam->getMinTimestep()) 
+                       && (dtdiff>(this->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 );
+                       rescale = true;
+                       if(!this->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);
+                       }
+               } // really change dt
+       }
+
+       if(mTimestepReduceLock>0) mTimestepReduceLock--;
+
+       
+       // forced back and forth switchting (for testing)
+       /*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) ){
+               rescale = true; minCutoff = false;
+               newdt = tadtogSwitch * this->mpParam->getTimestep();
+               this->mpParam->setDesiredTimestep( newdt );
+       } else 
+       if( ((this->mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) ||
+           ((this->mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){
+               rescale = true; minCutoff = false;
+               newdt = this->mpParam->getTimestep()/tadtogSwitch ;
+               this->mpParam->setDesiredTimestep( newdt );
+       } else {
+               rescale = false; minCutoff = false;
+       }
+       // */
+
+       // test mass rescale
+
+       scaleFac = newdt/this->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;
+
+               mTimeSwitchCounts++;
+               this->mpParam->calculateAllMissingValues( mSimulationTime, this->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();
+
+               // this might be called during init, before we have any particles
+               if(mpParticles) { mpParticles->adaptPartTimestep(scaleFac); }
+#if LBM_INCLUDE_TESTSOLVERS==1
+               if((mUseTestdata)&&(mpTest)) { 
+                       mpTest->adaptTimestep(scaleFac, mLevel[mMaxRefine].omega, mLevel[mMaxRefine].timestep, vec2L( mpParam->calculateGravity(mSimulationTime)) ); 
+                       mpTest->mGrav3d = mLevel[mMaxRefine].gravity;
+               }
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
+       
+               for(int lev=mMaxRefine; lev>=0 ; lev--) {
+                       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: "<<
+                                               " scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10);
+                       }
+                       if(lev!=mMaxRefine) coarseCalculateFluxareas(lev);
+
+                       int wss = 0, wse = 1;
+                       // only change currset (necessary for compressed grids, better for non-compr.gr.)
+                       wss = wse = mLevel[lev].setCurr;
+                       for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
+                                       // warning - check sets for higher levels...?
+                               FSGR_FORIJK1(lev) {
+                                       if( (RFLAG(lev,i,j,k, workSet) & CFBndMoving) ) {
+                                               /*
+                                               // paranoid check - shouldnt be necessary!
+                                               if(QCELL(lev, i, j, k, workSet, dFlux)!=mSimulationTime) {
+                                                       errMsg("TTTT","found invalid bnd cell.... removing at "<<PRINT_IJK);
+                                                       RFLAG(lev,i,j,k, workSet) = CFInter;
+                                                       // 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++) {
+                                                               QCELL(lev, i, j, k, workSet, l) = QCELL(lev, i, j, k, workSet, l)* scaleFac; 
+                                                       }
+                                               //} //  ok
+                                               continue;
+                                       }
+                                       if( 
+                                                       (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
+                                                       (RFLAG(lev,i,j,k, workSet) & CFInter) ||
+                                                       (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
+                                                       (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
+                                                       (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
+                                                       ) {
+                                               // these cells have to be scaled...
+                                       } else {
+                                               continue;
+                                       }
+
+                                       // collide on current set
+                                       LbmFloat rhoOld;
+                                       LbmVec velOld;
+                                       LbmFloat rho, ux,uy,uz;
+                                       rho=0.0; ux =  uy = uz = 0.0;
+                                       for(int l=0; l<this->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); 
+                                       } 
+                                       rhoOld = rho;
+                                       velOld = LbmVec(ux,uy,uz);
+
+                                       LbmFloat rhoNew = (rhoOld-rhoAvg)*scaleFac +rhoAvg;
+                                       LbmVec velNew = velOld * scaleFac;
+
+                                       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] );
+                                               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);
+                                       //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++) {
+                                               // org scaling
+                                               //df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important
+                                               // new scaling
+                                               LbmFloat dfn = feqNew[l] + (df[l]-feqOld[l])*dfScale*feqNew[l]/feqOld[l]; // non-eq. scaling, important
+                                               //df = eqNew + (df-eqOld)*dfScale; // modified ig scaling, no real difference?
+                                               QCELL(lev, i,j,k,workSet, l) = dfn;
+                                       }
+
+                                       if(RFLAG(lev,i,j,k, workSet) & CFInter) {
+                                               //if(workSet==mLevel[lev].setCurr) 
+                                               LbmFloat area = 1.0;
+                                               if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux);
+                                               massTOld += QCELL(lev, i,j,k,workSet, dMass) * area;
+                                               volTOld += QCELL(lev, i,j,k,workSet, dFfrac);
+
+                                               // wrong... QCELL(i,j,k,workSet, dMass] = (QCELL(i,j,k,workSet, dFfrac]*rhoNew);
+                                               QCELL(lev, i,j,k,workSet, dMass) = (QCELL(lev, i,j,k,workSet, dMass)/rhoOld*rhoNew);
+                                               QCELL(lev, i,j,k,workSet, dFfrac) = (QCELL(lev, i,j,k,workSet, dMass)/rhoNew);
+
+                                               //if(workSet==mLevel[lev].setCurr) 
+                                               massTNew += QCELL(lev, i,j,k,workSet, dMass);
+                                               volTNew += QCELL(lev, i,j,k,workSet, dFfrac);
+                                       }
+                                       if(RFLAG(lev,i,j,k, workSet) & CFFluid) { // DEBUG
+                                               if(RFLAG(lev,i,j,k, workSet) & (CFGrFromFine|CFGrFromCoarse)) { // DEBUG
+                                                       // dont include 
+                                               } else {
+                                                       LbmFloat area = 1.0;
+                                                       if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux) * mLevel[lev].lcellfactor;
+                                                       //if(workSet==mLevel[lev].setCurr) 
+                                                       massTOld += rhoOld*area;
+                                                       //if(workSet==mLevel[lev].setCurr) 
+                                                       massTNew += rhoNew*area;
+                                                       volTOld += area;
+                                                       volTNew += area;
+                                               }
+                                       }
+
+                               } // IJK
+                       } // workSet
+
+               } // lev
+
+               if(!this->mSilent) {
+                       debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<this->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);
+               }
+       } // 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()<<") " );
+               //brute force resacle all the time?
+
+               for(int lev=mMaxRefine; lev>=0 ; lev--) {
+               int rescs=0;
+               int wss = 0, wse = 1;
+#if COMPRESSGRIDS==1
+               if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
+#endif // COMPRESSGRIDS==1
+               for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
+               //for(int workSet = 0; workSet<=1; workSet++) {
+               FSGR_FORIJK1(lev) {
+
+                       //if( (RFLAG(lev, i,j,k, workSet) & CFFluid) || (RFLAG(lev, i,j,k, workSet) & CFInter) ) {
+                       if( 
+                                       (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
+                                       (RFLAG(lev,i,j,k, workSet) & CFInter) ||
+                                       (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
+                                       (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
+                                       (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
+                                       ) {
+                               // these cells have to be scaled...
+                       } else {
+                               continue;
+                       }
+
+                       // 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++) {
+                               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); 
+                       } 
+#ifndef WIN32
+                       if (!finite(rho)) {
+                               errMsg("adaptTimestep","Brute force non-finite rho at"<<PRINT_IJK);  // DEBUG!
+                               rho = 1.0;
+                               ux = uy = uz = 0.0;
+                               QCELL(lev, i, j, k, workSet, dMass) = 1.0;
+                               QCELL(lev, i, j, k, workSet, dFfrac) = 1.0;
+                       }
+#endif // WIN32
+
+                       if( (ux*ux+uy*uy+uz*uz)> (allowMax*allowMax) ) {
+                               LbmFloat cfac = allowMax/sqrt(ux*ux+uy*uy+uz*uz);
+                               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); }
+                               rescs++;
+                               debMsgDirect("B");
+                       }
+
+               } } 
+                       //if(rescs>0) { errMsg("adaptTimestep","!!!!! Brute force rescaling was necessary !!!!!!!"); }
+                       debMsgStd("adaptTimestep",DM_MSG,"Brute force rescale done. level:"<<lev<<" rescs:"<<rescs, 1);
+               //TTT mNumProblems += rescs; // add to problem display...
+               } // lev,set,ijk
+
+       } // try brute force rescale?
+
+       // time adap done...
+}
+
+
+
index b241834..10b172d 100644 (file)
@@ -208,8 +208,15 @@ class LbmFsgrSolver :
                virtual void parseAttrList();
                //! Initialize omegas and forces on all levels (for init/timestep change)
                void initLevelOmegas();
-               //! finish the init with config file values (allocate arrays...) 
-               virtual bool initializeSolver(); 
+
+               // multi step solver init
+               /*! finish the init with config file values (allocate arrays...) */
+               virtual bool initializeSolverMemory();
+               /*! init solver arrays */
+               virtual bool initializeSolverGrids();
+               /*! prepare actual simulation start, setup viz etc */
+               virtual bool initializeSolverPostinit();
+
                //! notify object that dump is in progress (e.g. for field dump) 
                virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
 
@@ -338,6 +345,8 @@ class LbmFsgrSolver :
                
                //! use time adaptivity? 
                bool mTimeAdap;
+               //! force smaller timestep for next LBM step? (eg for mov obj)
+               bool mForceTimeStepReduce;
 
                //! fluid vol height
                LbmFloat mFVHeight;
@@ -395,6 +404,12 @@ class LbmFsgrSolver :
                vector<LbmVec> mObjectSpeeds;
                //! partslip bc. values for obstacle boundary conditions
                vector<LbmFloat> mObjectPartslips;
+               //! moving object mass boundary condition values
+               vector<LbmFloat> mObjectMassMovnd;
+
+               //! permanent movobj vert storage
+         vector<ntlVec3Gfx>  mMOIVertices;
+       vector<ntlVec3Gfx>  mMOIVerticesOld;
 
                //! get isofield weights
                int mIsoWeightMethod;
@@ -459,9 +474,12 @@ class LbmFsgrSolver :
                void destroyTestdata();
                void handleTestdata();
                void set3dHeight(int ,int );
+
+               void initCpdata();
+               void handleCpdata();
        public:
-               // needed from testdata
-               void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy);
+               // needed for testdata
+               void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz);
 #endif // LBM_INCLUDE_TESTSOLVERS==1
 
        public: // former LbmModelLBGK  functions
@@ -473,7 +491,9 @@ class LbmFsgrSolver :
                inline void collideArrays( int i, int j, int k, // position - more for debugging
                                LbmFloat df[], LbmFloat &outrho, // out only!
                                // velocity modifiers (returns actual velocity!)
-                               LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, LbmFloat omega, LbmFloat csmago, LbmFloat *newOmegaRet, LbmFloat *newQoRet);
+                               LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
+                               LbmFloat omega, LbmVec gravity, LbmFloat csmago, 
+                               LbmFloat *newOmegaRet, LbmFloat *newQoRet);
 
 
                // former LBM models
index a417c33..908355f 100644 (file)
 // 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 )
+
+
 /*****************************************************************************/
 //! common variables 
 
@@ -296,7 +300,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mNumProblems(0), 
        mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
        mpPreviewSurface(NULL), 
-       mTimeAdap(true), 
+       mTimeAdap(true), mForceTimeStepReduce(false),
        mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
        mInitSurfaceSmoothing(0),
        mTimestepReduceLock(0),
@@ -305,11 +309,13 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mMinTimestep(0.0), mMaxTimestep(0.0),
        mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
        mDropMode(1), mDropSize(0.15), mDropSpeed(0.0),
-       mObjectSpeeds(), mObjectPartslips(),
+       mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
+       mMOIVertices(), mMOIVerticesOld(),
        mIsoWeightMethod(1),
        mMaxRefine(1), 
        mDfScaleUp(-1.0), mDfScaleDown(-1.0),
-       mInitialCsmago(0.04), mDebugOmegaRet(0.0),
+       mInitialCsmago(0.03), // set to 0.02 for mMaxRefine==0 below
+       mDebugOmegaRet(0.0),
        mLastOmega(1e10), mLastGravity(1e10),
        mNumInvIfTotal(0), mNumFsgrChanges(0),
        mDisableStandingFluidInit(0),
@@ -450,16 +456,6 @@ void LbmFsgrSolver::parseAttrList()
        this->mSmoothSurface = this->mpAttrs->readFloat("smoothsurface", this->mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
        this->mSmoothNormals = this->mpAttrs->readFloat("smoothnormals", this->mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
 
-       mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
-       // deprecated!
-       float mInitialCsmagoCoarse = 0.0;
-       mInitialCsmagoCoarse = this->mpAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
-#if USE_LES==1
-#else // USE_LES==1
-       debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
-       mInitialCsmago = 0.0;
-#endif // USE_LES==1
-
        // refinement
        mMaxRefine = this->mRefinementDesired;
        mMaxRefine  = this->mpAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
@@ -481,12 +477,24 @@ void LbmFsgrSolver::parseAttrList()
        mUseTestdata=1; // DEBUG
 #endif // ELBEEM_PLUGIN
        errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
+
 #else // LBM_INCLUDE_TESTSOLVERS!=1
        // off by default
        mUseTestdata = 0;
        if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
-       if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
 #endif // LBM_INCLUDE_TESTSOLVERS!=1
+       if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
+
+       if(mMaxRefine==0) mInitialCsmago=0.02;
+       mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
+       // deprecated!
+       float mInitialCsmagoCoarse = 0.0;
+       mInitialCsmagoCoarse = this->mpAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
+#if USE_LES==1
+#else // USE_LES==1
+       debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
+       mInitialCsmago = 0.0;
+#endif // USE_LES==1
 }
 
 
@@ -571,7 +579,9 @@ void LbmFsgrSolver::initLevelOmegas()
 /******************************************************************************
  * Init Solver (values should be read from config file)
  *****************************************************************************/
-bool LbmFsgrSolver::initializeSolver()
+
+/*! finish the init with config file values (allocate arrays...) */
+bool LbmFsgrSolver::initializeSolverMemory()
 {
   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<this->mInitDone<<" "<<(void*)this,1);
 
@@ -597,7 +607,7 @@ bool LbmFsgrSolver::initializeSolver()
                calculateMemreqEstimate( this->mSizex, this->mSizey, this->mSizez, mMaxRefine, &memEstFromFunc, &memreqStr );
                
                double memLimit;
-               if(sizeof(int)==4) {
+               if(sizeof(void*)==4) {
                        // 32bit system, limit to 2GB
                        memLimit = 2.0* 1024.0*1024.0*1024.0;
                } else {
@@ -619,6 +629,7 @@ bool LbmFsgrSolver::initializeSolver()
        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*) <<
          ", est. Mem.Req.: "<<memreqStr        ,2);
        this->mpParam->setSize(this->mSizex, this->mSizey, this->mSizez);
 
@@ -835,10 +846,43 @@ bool LbmFsgrSolver::initializeSolver()
                }
        }
 
+       if(LBMDIM==2) {
+               if(this->mOutputSurfacePreview) {
+                       errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
+                       this->mOutputSurfacePreview = 0; }
+       }
+#if LBM_USE_GUI==1
+       if(this->mOutputSurfacePreview) {
+               errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
+               this->mOutputSurfacePreview = 0; }
+#endif // LBM_USE_GUI==1
+       if(this->mOutputSurfacePreview) {
+               // same as normal one, but use reduced size
+               mpPreviewSurface = new IsoSurface( this->mIsoValue );
+               mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
+               mpPreviewSurface->setIsolevel( this->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" );
+               mpPreviewSurface->setName( "preview" );
+       
+               debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*this->mSizex)<<","<<(pfac*this->mSizey)<<","<<(pfac*this->mSizez)<<" enabled",10);
+       }
+
        // init defaults
        mAvgNumUsedCells = 0;
        this->mFixMass= 0.0;
+       return true;
+}
 
+/*! init solver arrays */
+bool LbmFsgrSolver::initializeSolverGrids() {
   /* init boundaries */
   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
        // init obstacles, and reinit time step size 
@@ -919,6 +963,8 @@ bool LbmFsgrSolver::initializeSolver()
        }
        // Symmetry tests */
 
+       //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
+
        // prepare interface cells
        initFreeSurfaces();
        initStandingFluidGradient();
@@ -956,9 +1002,12 @@ bool LbmFsgrSolver::initializeSolver()
                }
        }
 #endif // ELBEEM_PLUGIN!=1
-       
+       return true;
+}
 
-       // ----------------------------------------------------------------------
+
+/*! prepare actual simulation start, setup viz etc */
+bool LbmFsgrSolver::initializeSolverPostinit() {
        // coarsen region
        myTime_t fsgrtstart = getTime(); 
        for(int lev=mMaxRefine-1; lev>=0; lev--) {
@@ -991,38 +1040,8 @@ bool LbmFsgrSolver::initializeSolver()
        } } // first copy flags */
 
 
-       
-       if(LBMDIM==2) {
-               if(this->mOutputSurfacePreview) {
-                       errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
-                       this->mOutputSurfacePreview = 0; }
-       }
-#if LBM_USE_GUI==1
-       if(this->mOutputSurfacePreview) {
-               errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
-               this->mOutputSurfacePreview = 0; }
-#endif // LBM_USE_GUI==1
-       
-       if(this->mOutputSurfacePreview) {
-
-               // same as normal one, but use reduced size
-               mpPreviewSurface = new IsoSurface( this->mIsoValue );
-               mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
-               mpPreviewSurface->setIsolevel( this->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" );
-               mpPreviewSurface->setName( "preview" );
-       
-               debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*this->mSizex)<<","<<(pfac*this->mSizey)<<","<<(pfac*this->mSizez)<<" enabled",10);
-       }
-
+       // old mpPreviewSurface init
+       //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
        // make sure fill fracs are right for first surface generation
        stepMain();
 
@@ -1036,11 +1055,12 @@ bool LbmFsgrSolver::initializeSolver()
 
 
        // now really done...
-  debugOut("LbmFsgrSolver::initialize : SurfaceGen: SmsOrg("<<this->mSmoothSurface<<","<<this->mSmoothNormals<<","<<featureSize<<"), Iso("<<this->mpIso->getSmoothSurface()<<","<<this->mpIso->getSmoothNormals()<<") ",10);
+  debugOut("LbmFsgrSolver::initialize : SurfaceGen: SmsOrg("<<this->mSmoothSurface<<","<<this->mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<this->mpIso->getSmoothSurface()<<","<<this->mpIso->getSmoothNormals()<<") ",10);
   debugOut("LbmFsgrSolver::initialize : Init done ... ",10);
        this->mInitDone = 1;
 
 #if LBM_INCLUDE_TESTSOLVERS==1
+       initCpdata();
        initTestdata();
 #endif // ELBEEM_PLUGIN!=1
        // not inited? dont use...
@@ -1052,14 +1072,54 @@ bool LbmFsgrSolver::initializeSolver()
 
 
 
+// macros for mov obj init
+#if LBMDIM==2
+
+#define POS2GRID_CHECK(vec,n) \
+                               monTotal++;\
+                               const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
+                               if(i<=0) continue; \
+                               if(i>=mLevel[level].lSizex-1) continue; \
+                               const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
+                               if(j<=0) continue; \
+                               if(j>=mLevel[level].lSizey-1) continue;  \
+                               const int k=0; \
+
+#else // LBMDIM -> 3
+#define POS2GRID_CHECK(vec,n) \
+                               monTotal++;\
+                               const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
+                               if(i<=0) continue; \
+                               if(i>=mLevel[level].lSizex-1) continue; \
+                               const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
+                               if(j<=0) continue; \
+                               if(j>=mLevel[level].lSizey-1) continue; \
+                               const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
+                               if(k<=0) continue; \
+                               if(k>=mLevel[level].lSizez-1) continue; \
+
+#endif // LBMDIM 
+
+// calculate object velocity from vert arrays in objvel vec
+#define OBJVEL_CALC \
+                               LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
+                               const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
+                               USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
+                               if(usqr>maxusqr) { \
+                                       /* cutoff at maxVelVal */ \
+                                       for(int jj=0; jj<3; jj++) { \
+                                               if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
+                                               if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
+                                       } \
+                               } }
 
 /*****************************************************************************/
 //! init moving obstacles for next sim step sim 
 /*****************************************************************************/
 void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
-
        myTime_t monstart = getTime();
-       // new test
+
+       // movobj init
        const int level = mMaxRefine;
        const int workSet = mLevel[level].setCurr;
        LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
@@ -1072,6 +1132,9 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
        //if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);
 
+       const LbmFloat maxVelVal = 0.1666;
+       const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
+
        LbmFloat rhomass = 0.0;
        CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
        CellFlagType ntype = CFInvalid;
@@ -1088,30 +1151,53 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                //iniPos[2] = this->mvGeoStart[2] + dvec[2]*getForZMin1();
        }
        
+       if( (int)mObjectMassMovnd.size() < numobjs) {
+               for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
+                       mObjectMassMovnd.push_back(0.);
+               }
+       }
+       
        // stats
        int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
-       CellFlagType nbflag[LBM_DFNUM]; 
        int nbored;
-  vector<ntlVec3Gfx>  vertices;
-  vector<ntlVec3Gfx>  verticesOld;
-       int i,j,k; // grid pos init
-       LbmFloat ux,uy,uz, rho;
        for(int OId=0; OId<numobjs; OId++) {
                ntlGeometryObject *obj = (*this->mpGiObjects)[OId];
-               if( (!staticInit) && (!obj->getIsAnimated()) ) continue;
-               if( ( staticInit) && ( obj->getIsAnimated()) ) continue;
+               bool skip = false;
+               if(obj->getGeoInitId() != this->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;
 
-               if(obj->getGeoInitType()&FGI_ALLBOUNDS) {
+               if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || 
+                               (obj->getGeoInitType()&FGI_FLUID) && staticInit ) {
 
                        otype = ntype = CFInvalid;
                        switch(obj->getGeoInitType()) {
                                case FGI_BNDPART: 
                                case FGI_BNDFREE: 
-                                       errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
-                               case FGI_BNDNO: 
-                                       rhomass = BND_FILL;
+                                       if(!staticInit) {
+                                               errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
+                                               otype = ntype = CFBnd|CFBndNoslip;
+                                       } else {
+                                               if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
+                                               if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
+                                       }
+                                       break; 
+                                       /*
+                               case FGI_BNDPART: rhomass = BND_FILL;
+                                       otype = ntype = CFBnd|CFBndPartslip;
+                                       break;
+                               case FGI_BNDFREE: rhomass = BND_FILL;
+                                       otype = ntype = CFBnd|CFBndFreeslip;
+                                       break;
+                               // */
+                               case FGI_BNDNO:   rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndNoslip;
                                        break;
+                               case FGI_FLUID: 
+                                       otype = ntype = CFFluid; 
+                                       break;
                                case FGI_MBNDINFLOW: 
                                        otype = ntype = CFMbndInflow; 
                                        break;
@@ -1123,23 +1209,35 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                        int active =    ((obj->getGeoActive(targetTime)>0.)? 1:0);
                        //errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<"  s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
                        // skip inactive in/out flows
+                       if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
                        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 );
 
-                       vertices.clear();
-                       obj->getMovingPoints(vertices);
-                       verticesOld = vertices;
-                       // WARNING - assumes mSimulationTime is global!?
-                       obj->applyTransformation(targetTime, &vertices,NULL, 0, vertices.size(), false );
-                       monTrafo += vertices.size();
-
-                       // correct flags from last position, but extrapolate
-                       // velocity to next timestep
-                       obj->applyTransformation(sourceTime, &verticesOld,NULL, 0, verticesOld.size(), false );
-                       monTrafo += verticesOld.size();
+                       mMOIVertices.clear();
+                       if(obj->getMeshAnimated()) { 
+                               // do two full update
+                               mMOIVerticesOld.clear();
+                               obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices,  mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
+                               monTrafo += mMOIVerticesOld.size();
+                               obj->applyTransformation(sourceTime, &mMOIVerticesOld,NULL, 0, mMOIVerticesOld.size(), false );
+                               monTrafo += mMOIVertices.size();
+                               obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
+                       } else {
+                               // only do transform update
+                               obj->getMovingPoints(mMOIVertices);
+                               mMOIVerticesOld = mMOIVertices;
+                               // WARNING - assumes mSimulationTime is global!?
+                               obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
+                               monTrafo += mMOIVertices.size();
+
+                               // correct flags from last position, but extrapolate
+                               // velocity to next timestep
+                               obj->applyTransformation(sourceTime, &mMOIVerticesOld,NULL, 0, mMOIVerticesOld.size(), false );
+                               monTrafo += mMOIVerticesOld.size();
+                       }
 
                        // object types
                        if(ntype&CFBnd){
@@ -1153,144 +1251,189 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                        ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - this->mpParam->getTimestep(),sourceTime);
                                        if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
                                }
-
-#if LBMDIM==2
-#define CHECKIJK \
-                               if(i<=0) continue; \
-                               if(j<=0) continue; \
-                               if(i>=mLevel[level].lSizex-2) continue; \
-                               if(j>=mLevel[level].lSizey-2) continue; 
-#define POS2GRID(vec,n) \
-                               monTotal++;\
-                               i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
-                               j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
-                               k=0;
-#else // LBMDIM -> 3
-#define CHECKIJK \
-                               if(i<=0) continue; \
-                               if(j<=0) continue; \
-                               if(k<=0) continue; \
-                               if(i>=mLevel[level].lSizex-2) continue; \
-                               if(j>=mLevel[level].lSizey-2) continue; \
-                               if(k>=mLevel[level].lSizez-2) continue; 
-#define POS2GRID(vec,n) \
-                               monTotal++;\
-                               i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
-                               j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
-                               k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); 
-#endif // LBMDIM 
+                               if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
+                               CellFlagType rflagnb[27];
+                               LbmFloat massCheck = 0.;
+                               int massReinits=0;                              
+                               bool fillCells = (mObjectMassMovnd[OId]<=-1.);
+                               LbmFloat massCScale; //, massCScalePos, massCScaleNeg;
+                               //if(mInitialMass>0.) {
+                                       //massCScale = (mInitialMass-mObjectMassMovnd[OId])/mInitialMass;
+                                       massCScale = 1.-(mObjectMassMovnd[OId]/(LbmFloat)(mMOIVertices.size()/10) );