Only restructured code:

- added Hos fixes
- split up solver into 3 cpp files (as suggested by jonathan)
- de-inlined function that caused gcc33 to use >1GB of memory

intern/elbeem/SConscript

@@ -29,22 +29,25 @@ else:
 
 		"intern/attributes.cpp",
 		"intern/elbeem.cpp",
-		"intern/factory_fsgr.cpp",
 		"intern/isosurface.cpp",
-		"intern/lbminterface.cpp",
 		"intern/ntl_blenderdumper.cpp",
 		"intern/ntl_bsptree.cpp",
 		"intern/ntl_geometrymodel.cpp",
 		"intern/ntl_geometryobject.cpp",
 		"intern/ntl_lightobject.cpp",
 		"intern/ntl_ray.cpp",
-		"intern/ntl_raytracer.cpp",
 		"intern/ntl_scene.cpp",
+		"intern/ntl_world.cpp",
 		"intern/parametrizer.cpp",
 		"intern/particletracer.cpp",
 		"intern/simulation_object.cpp",
 		"intern/utilities.cpp",
-		"intern/blendercall.cpp"
+		"intern/blendercall.cpp",
+
+		"intern/solver_init.cpp",
+		"intern/solver_interface.cpp",
+		"intern/solver_main.cpp",
+		"intern/solver_util.cpp"
 
 		]; # sources
 

intern/elbeem/intern/blendercall.cpp

@@ -10,7 +10,6 @@
  *****************************************************************************/
  
 #include "globals.h"
-#include "ntl_raytracer.h"
 #include "ntl_blenderdumper.h"
 #include <stdlib.h>
 

intern/elbeem/intern/cfglexer.cpp

@@ -810,9 +810,9 @@ char *yy_text;
 #define CHAR_BUFFER_SIZE 8000
 char charBuffer[ CHAR_BUFFER_SIZE ];
 
-extern "C" int yy_wrap (void ) { return 1; }
 int lineCount = 1;
 
+extern "C" int yy_wrap (void ) { return 1; }
 #define YY_NO_UNISTD_H
 	
 /*----------------------------------------------------------------------------*/

intern/elbeem/intern/cfgparser.cpp

@@ -1233,7 +1233,6 @@ yysymprint (yyoutput, yytype, yyvaluep)
 # endif
   switch (yytype)
     {
-			case 0:
       default:
         break;
     }

intern/elbeem/intern/elbeem.cpp

@@ -21,8 +21,8 @@ double guiRoiEZ = 1.0;
 int guiRoiMaxLev=6, guiRoiMinLev=0;
 
 //! global raytracer pointer (=world)
-class ntlRaytracer;
+class ntlWorld;
-ntlRaytracer *gpWorld = (ntlRaytracer*)0;
+ntlWorld *gpWorld = (ntlWorld*)0;
 
 //! debug output switch
 bool myDebugOut = false;

intern/elbeem/intern/factory_fsgr.cpp

@@ -1,29 +0,0 @@
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * Copyright 2003,2004,2005 Nils Thuerey
- *
- * Standard LBM Factory implementation
- *
- *****************************************************************************/
-
-// disable sometimes to speed up compiling/2d tests
-#define DISABLE 0
-
-#include "lbmfsgrsolver.h"
-#include "factory_lbm.h"
-
-//! lbm factory functions
-LbmSolverInterface* createSolverLbmFsgr() {
-#if DISABLE!=1
-#if LBMDIM==2
-	return new LbmFsgrSolver< LbmBGK2D >();
-#endif // LBMDIM==2
-#if LBMDIM==3
-	return new LbmFsgrSolver< LbmBGK3D >();
-#endif // LBMDIM==3
-#endif // DISABLE
-	return NULL;
-}
-
-

intern/elbeem/intern/factory_lbm.h

@@ -1,18 +0,0 @@
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * Copyright 2003,2004 Nils Thuerey
- *
- * 2D/3D LBM Factory header
- *
- *****************************************************************************/
-
-#include "lbminterface.h"
-
-//! lbm factory functions
-LbmSolverInterface* createSolverLbmFsgr();
-#ifdef LBM_INCLUDE_TESTSOLVERS
-LbmSolverInterface* createSolverOld();
-#endif // LBM_INCLUDE_OLD
-
-

intern/elbeem/intern/globals.h

@@ -15,8 +15,8 @@
 //! user interface variables
 
 // global raytracer pointer (=world)
-class ntlRaytracer;
+class ntlWorld;
-extern ntlRaytracer *gpWorld;
+extern ntlWorld *gpWorld;
 
 // debug output switch
 extern bool myDebugOut;

intern/elbeem/intern/isosurface.cpp

@@ -14,10 +14,13 @@
 #include <algorithm>
 #include <stdio.h>
 
+
+// sirdude fix for solaris
 #if !defined(linux) && (defined (__sparc) || defined (__sparc__))
 #include <ieeefp.h>
 #endif
 
+
 /******************************************************************************
  * Constructor
  *****************************************************************************/
@@ -176,9 +179,6 @@ void IsoSurface::triangulate( void )
 				value[5] = *getData(i+1,j  ,k+1);
 				value[6] = *getData(i+1,j+1,k+1);
 				value[7] = *getData(i  ,j+1,k+1);
-				//errMsg("ISOT2D"," at "<<PRINT_IJK<<" "
-						//<<" v0="<<value[0] <<" v1="<<value[1] <<" v2="<<value[2] <<" v3="<<value[3]
-						//<<" v4="<<value[4] <<" v5="<<value[5] <<" v6="<<value[6] <<" v7="<<value[7] );
 
 				// check intersections of isosurface with edges, and calculate cubie index
 				cubeIndex = 0;
@@ -235,34 +235,7 @@ void IsoSurface::triangulate( void )
 							const ntlVec3Gfx p2 = pos[ e2  ];    // scalar field pos 2
 							const float valp1  = value[ e1  ];  // scalar field val 1
 							const float valp2  = value[ e2  ];  // scalar field val 2
-							//double mu;            // interpolation value
-							//ntlVec3Gfx p;         // new point
 
-							// choose if point should be calculated by interpolation,
-							// or "Karolin" method 
-
-							//double deltaVal = ABS(valp2-valp1);
-							//if(deltaVal >-10.0) {
-							// standard interpolation
-							//vertInfo[e].type = 0; 
-							/*if (ABS(mIsoValue-valp1) < 0.000000001) {
-								mu = 0.0;
-							} else {
-								if (ABS(mIsoValue-valp2) < 0.000000001) {
-									mu = 1.0;
-								} else {
-									mu = (mIsoValue - valp1) / (valp2 - valp1);
-								}
-							} */
-							/*} else {
-								errorOut(" ? ");
-							// use fill grade (=karo) 
-							vertInfo[e].type = 1; 
-							if(valp1 < valp2) { mu = 1.0- (valp1 + valp2 - mIsoValue);
-							} else { mu = 0.0+ (valp1 + valp2 - mIsoValue); } 
-							} */
-
-							//const float mu = (mIsoValue - valp1) / (valp2 - valp1);
 							float mu;
 							if(valp1 < valp2) {
 								mu = 1.0 - 1.0*(valp1 + valp2 - mIsoValue);
@@ -270,15 +243,11 @@ void IsoSurface::triangulate( void )
 								mu = 0.0 + 1.0*(valp1 + valp2 - mIsoValue);
 							}
 
-							float isov2 = mIsoValue;
+							//float isov2 = mIsoValue;
-							isov2 = (valp1+valp2)*0.5;
+							//isov2 = (valp1+valp2)*0.5;
-							mu = (isov2 - valp1) / (valp2 - valp1);
+							//mu = (isov2 - valp1) / (valp2 - valp1);
-							mu = (isov2) / (valp2 - valp1);
+							//mu = (isov2) / (valp2 - valp1);
-
 							mu = (mIsoValue - valp1) / (valp2 - valp1);
-							//mu *= mu;
-
-
 
 							// init isolevel vertex
 							ilv.v = p1 + (p2-p1)*mu;
@@ -319,14 +288,9 @@ void IsoSurface::triangulate( void )
 		normalize( mPoints[ni].n );
 	}
 
-	if(mSmoothSurface>0.0) {
+	//for(int i=0; i<mLoopSubdivs; i++) { subdivide(); }// DEBUG test
-		// not needed for post normal smoothing? 
+	if(mSmoothSurface>0.0) { smoothSurface(mSmoothSurface); }
-		// if(mSmoothNormals<=0.0) { smoothNormals(mSmoothSurface*0.5); }
+	if(mSmoothNormals>0.0) { smoothNormals(mSmoothNormals); }
-		smoothSurface(mSmoothSurface);
-	}
-	if(mSmoothNormals>0.0) {
-		smoothNormals(mSmoothNormals);
-	}
 	myTime_t tritimeend = getTime(); 
 	debMsgStd("IsoSurface::triangulate",DM_MSG,"took "<< ((tritimeend-tritimestart)/(double)1000.0)<<"s, ss="<<mSmoothSurface<<" sm="<<mSmoothNormals , 10 );
 }
@@ -444,6 +408,283 @@ inline ntlVec3Gfx IsoSurface::getNormal(int i, int j,int k) {
  * 
  *****************************************************************************/
 
+
+// Subdivide a mesh allways loop
+/*void IsoSurface::subdivide()
+{
+	int i;
+
+	mAdjacentFaces.clear(); 
+	mAcrossEdge.clear();
+
+	//void TriMesh::need_adjacentfaces()
+	{
+		vector<int> numadjacentfaces(mPoints.size());
+		//errMsg("SUBDIV ADJFA1", " "<<mPoints.size()<<" - "<<numadjacentfaces.size() );
+		int i;
+		for (i = 0; i < (int)mIndices.size()/3; i++) {
+			numadjacentfaces[mIndices[i*3 + 0]]++;
+			numadjacentfaces[mIndices[i*3 + 1]]++;
+			numadjacentfaces[mIndices[i*3 + 2]]++;
+		}
+
+		mAdjacentFaces.resize(mPoints.size());
+		for (i = 0; i < (int)mPoints.size(); i++)
+			mAdjacentFaces[i].reserve(numadjacentfaces[i]);
+
+		for (i = 0; i < (int)mIndices.size()/3; i++) {
+			for (int j = 0; j < 3; j++)
+				mAdjacentFaces[mIndices[i*3 + j]].push_back(i);
+		}
+
+	}
+
+	// Find the face across each edge from each other face (-1 on boundary)
+	// If topology is bad, not necessarily what one would expect...
+	//void TriMesh::need_across_edge()
+	{
+		mAcrossEdge.resize(mIndices.size(), -1);
+
+		for (int i = 0; i < (int)mIndices.size()/3; i++) {
+			for (int j = 0; j < 3; j++) {
+				if (mAcrossEdge[i*3 + j] != -1)
+					continue;
+				int v1 = mIndices[i*3 + ((j+1)%3)];
+				int v2 = mIndices[i*3 + ((j+2)%3)];
+				const vector<int> &a1 = mAdjacentFaces[v1];
+				const vector<int> &a2 = mAdjacentFaces[v2];
+				for (int k1 = 0; k1 < (int)a1.size(); k1++) {
+					int other = a1[k1];
+					if (other == i)
+						continue;
+					vector<int>::const_iterator it =
+						std::find(a2.begin(), a2.end(), other);
+					if (it == a2.end())
+						continue;
+
+					//int ind = (faces[other].indexof(v1)+1)%3;
+					int ind = -1;
+					if( mIndices[other*3+0] == (unsigned int)v1 ) ind = 0;
+					else if( mIndices[other*3+1] == (unsigned int)v1 ) ind = 1;
+					else if( mIndices[other*3+2] == (unsigned int)v1 ) ind = 2;
+					ind = (ind+1)%3;
+
+					if ( (int)mIndices[other*3 + ((ind+1)%3)] != v2)
+						continue;
+					mAcrossEdge[i*3 + j] = other;
+					mAcrossEdge[other*3 + ind] = i;
+					break;
+				}
+			}
+		}
+
+		//errMsg("SUBDIV ACREDG", "Done.\n");
+	}
+
+	//errMsg("SUBDIV","start");
+	// Introduce new vertices
+	int nf = (int)mIndices.size() / 3;
+
+	//vector<TriMesh::Face> newverts(nf, TriMesh::Face(-1,-1,-1));
+	vector<int> newverts(nf*3); //, TriMesh::Face(-1,-1,-1));
+	for(int j=0; j<(int)newverts.size(); j++) newverts[j] = -1;
+
+	int old_nv = (int)mPoints.size();
+	mPoints.reserve(4 * old_nv);
+	vector<int> newvert_count(old_nv + 3*nf); // wichtig...?
+	//errMsg("NC", newvert_count.size() );
+
+	for (i = 0; i < nf; i++) {
+		for (int j = 0; j < 3; j++) {
+			int ae = mAcrossEdge[i*3 + j];
+			if (newverts[i*3 + j] == -1 && ae != -1) {
+				if (mAcrossEdge[ae*3 + 0] == i)
+					newverts[i*3 + j] = newverts[ae*3 + 0];
+				else if (mAcrossEdge[ae*3 + 1] == i)
+					newverts[i*3 + j] = newverts[ae*3 + 1];
+				else if (mAcrossEdge[ae*3 + 2] == i)
+					newverts[i*3 + j] = newverts[ae*3 + 2];
+			}
+			if (newverts[i*3 + j] == -1) {
+				IsoLevelVertex ilv;
+				ilv.v = ntlVec3Gfx(0.0);
+				ilv.n = ntlVec3Gfx(0.0);
+				mPoints.push_back(ilv);
+				newverts[i*3 + j] = (int)mPoints.size() - 1;
+				if (ae != -1) {
+					if (mAcrossEdge[ae*3 + 0] == i)
+						newverts[ae*3 + 0] = newverts[i*3 + j];
+					else if (mAcrossEdge[ae*3 + 1] == i)
+						newverts[ae*3 + 1] = newverts[i*3 + j];
+					else if (mAcrossEdge[ae*3 + 2] == i)
+						newverts[ae*3 + 2] = newverts[i*3 + j];
+				}
+			}
+			if(ae != -1) {
+				mPoints[newverts[i*3 + j]].v +=
+					mPoints[ mIndices[i*3 + ( j     )] ].v * 0.25f  + // j = 0,1,2?
+					mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.375f +
+					mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.375f;
+#if RECALCNORMALS==0
+				mPoints[newverts[i*3 + j]].n +=
+					mPoints[ mIndices[i*3 + ( j     )] ].n * 0.25f  + // j = 0,1,2?
+					mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.375f +
+					mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.375f;
+#endif // RECALCNORMALS==0
+			} else {
+				mPoints[newverts[i*3 + j]].v +=
+					mPoints[ mIndices[i*3 + ((j+1)%3)] ].v * 0.5f   +
+					mPoints[ mIndices[i*3 + ((j+2)%3)] ].v * 0.5f  ;
+#if RECALCNORMALS==0
+				mPoints[newverts[i*3 + j]].n +=
+					mPoints[ mIndices[i*3 + ((j+1)%3)] ].n * 0.5f   +
+					mPoints[ mIndices[i*3 + ((j+2)%3)] ].n * 0.5f  ;
+#endif // RECALCNORMALS==0
+			}
+
+			newvert_count[newverts[i*3 + j]]++;
+		}
+	}
+	for (i = old_nv; i < (int)mPoints.size(); i++) {
+		if (!newvert_count[i])
+			continue;
+		float scale = 1.0f / newvert_count[i];
+		mPoints[i].v *= scale;
+
+#if RECALCNORMALS==0
+		//mPoints[i].n *= scale;
+		//normalize( mPoints[i].n );
+#endif // RECALCNORMALS==0
+	}
+
+	// Update old vertices
+	for (i = 0; i < old_nv; i++) {
+			ntlVec3Gfx bdyavg(0.0), nbdyavg(0.0);
+			ntlVec3Gfx norm_bdyavg(0.0), norm_nbdyavg(0.0); // N
+			int nbdy = 0, nnbdy = 0;
+			int naf = (int)mAdjacentFaces[i].size();
+			if (!naf)
+				continue;
+			for (int j = 0; j < naf; j++) {
+				int af = mAdjacentFaces[i][j];
+
+				int afi = -1;
+				if( mIndices[af*3+0] == (unsigned int)i ) afi = 0;
+				else if( mIndices[af*3+1] == (unsigned int)i ) afi = 1;
+				else if( mIndices[af*3+2] == (unsigned int)i ) afi = 2;
+
+				int n1 = (afi+1) % 3;
+				int n2 = (afi+2) % 3;
+				if (mAcrossEdge[af*3 + n1] == -1) {
+					bdyavg += mPoints[newverts[af*3 + n1]].v;
+#if RECALCNORMALS==0
+					//norm_bdyavg += mPoints[newverts[af*3 + n1]].n;
+#endif // RECALCNORMALS==0
+					nbdy++;
+				} else {
+					nbdyavg += mPoints[newverts[af*3 + n1]].v;
+#if RECALCNORMALS==0
+					//norm_nbdyavg += mPoints[newverts[af*3 + n1]].n;
+#endif // RECALCNORMALS==0
+					nnbdy++;
+				}
+				if (mAcrossEdge[af*3 + n2] == -1) {
+					bdyavg += mPoints[newverts[af*3 + n2]].v;
+#if RECALCNORMALS==0
+					//norm_bdyavg += mPoints[newverts[af*3 + n2]].n;
+#endif // RECALCNORMALS==0
+					nbdy++;
+				} else {
+					nbdyavg += mPoints[newverts[af*3 + n2]].v;
+#if RECALCNORMALS==0
+					//norm_nbdyavg += mPoints[newverts[af*3 + n2]].n;
+#endif // RECALCNORMALS==0
+					nnbdy++;
+				}
+			}
+
+			float alpha;
+			ntlVec3Gfx newpt;
+			if (nbdy) {
+				newpt = bdyavg / (float) nbdy;
+				alpha = 0.5f;
+			} else if (nnbdy) {
+				newpt = nbdyavg / (float) nnbdy;
+				if (nnbdy == 6)
+					alpha = 1.05;
+				else if (nnbdy == 8)
+					alpha = 0.86;
+				else if (nnbdy == 10)
+					alpha = 0.7;
+				else
+					alpha = 0.6;
+			} else {
+				continue;
+			}
+			mPoints[i].v *= 1.0f - alpha;
+			mPoints[i].v += newpt * alpha;
+
+#if RECALCNORMALS==0
+			//mPoints[i].n *= 1.0f - alpha;
+			//mPoints[i].n += newpt * alpha;
+#endif // RECALCNORMALS==0
+	}
+
+	// Insert new faces
+	mIndices.reserve(4*nf);
+	for (i = 0; i < nf; i++) {
+		mIndices.push_back( mIndices[i*3 + 0]);
+		mIndices.push_back( newverts[i*3 + 2]);
+		mIndices.push_back( newverts[i*3 + 1]);
+
+		mIndices.push_back( mIndices[i*3 + 1]);
+		mIndices.push_back( newverts[i*3 + 0]);
+		mIndices.push_back( newverts[i*3 + 2]);
+
+		mIndices.push_back( mIndices[i*3 + 2]);
+		mIndices.push_back( newverts[i*3 + 1]);
+		mIndices.push_back( newverts[i*3 + 0]);
+
+		mIndices[i*3+0] = newverts[i*3+0];
+		mIndices[i*3+1] = newverts[i*3+1];
+		mIndices[i*3+2] = newverts[i*3+2];
+	}
+
+	// recalc normals
+#if RECALCNORMALS==1
+	{
+		int nf = (int)mIndices.size()/3, nv = (int)mPoints.size(), i;
+		for (i = 0; i < nv; i++) {
+			mPoints[i].n = ntlVec3Gfx(0.0);
+		}
+		for (i = 0; i < nf; i++) {
+			const ntlVec3Gfx &p0 = mPoints[mIndices[i*3+0]].v;
+			const ntlVec3Gfx &p1 = mPoints[mIndices[i*3+1]].v;
+			const ntlVec3Gfx &p2 = mPoints[mIndices[i*3+2]].v;
+			ntlVec3Gfx a = p0-p1, b = p1-p2, c = p2-p0;
+			float l2a = normNoSqrt(a), l2b = normNoSqrt(b), l2c = normNoSqrt(c);
+
+			ntlVec3Gfx facenormal = cross(a, b);
+
+			mPoints[mIndices[i*3+0]].n += facenormal * (1.0f / (l2a * l2c));
+			mPoints[mIndices[i*3+1]].n += facenormal * (1.0f / (l2b * l2a));
+			mPoints[mIndices[i*3+2]].n += facenormal * (1.0f / (l2c * l2b));
+		}
+
+		for (i = 0; i < nv; i++) {
+			normalize(mPoints[i].n);
+		}
+	}
+#else // RECALCNORMALS==1
+		for (i = 0; i < (int)mPoints.size(); i++) {
+			normalize(mPoints[i].n);
+		}
+#endif // RECALCNORMALS==1
+
+	//errMsg("SUBDIV","done nv:"<<mPoints.size()<<" nf:"<<mIndices.size() );
+}*/
+
+
 // Diffuse a vector field at 1 vertex, weighted by
 // a gaussian of width 1/sqrt(invsigma2)
 void IsoSurface::diffuseVertexField(ntlVec3Gfx *field, const int pointerScale, int v, float invsigma2, ntlVec3Gfx &flt)

intern/elbeem/intern/lbmdimensions.h

@@ -1,391 +0,0 @@
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * All code distributed as part of El'Beem is covered by the version 2 of the 
- * GNU General Public License. See the file COPYING for details.
- * Copyright 2003-2005 Nils Thuerey
- *
- * Combined 2D/3D Lattice Boltzmann Solver auxiliary classes
- * 
- *****************************************************************************/
-#ifndef LBMHEADER_H
-
-/* LBM Files */
-#include "lbminterface.h"
-#include <sstream>
-
-
-//! shorten static const definitions
-#define STCON static const
-
-
-/*****************************************************************************/
-/*! class for solver templating - 3D implementation */
-//class LbmD3Q19 : public LbmSolverInterface {
-class LbmD3Q19 {
-
-	public:
-
-		// constructor, init interface
-		LbmD3Q19() {};
-		// virtual destructor 
-		virtual ~LbmD3Q19() {};
-		//! id string of solver
-		string getIdString() { return string("3D"); }
-
-		//! how many dimensions?
-		STCON int cDimension;
-
-		// Wi factors for collide step 
-		STCON LbmFloat cCollenZero;
-		STCON LbmFloat cCollenOne;
-		STCON LbmFloat cCollenSqrtTwo;
-
-		//! threshold value for filled/emptied cells 
-		STCON LbmFloat cMagicNr2;
-		STCON LbmFloat cMagicNr2Neg;
-		STCON LbmFloat cMagicNr;
-		STCON LbmFloat cMagicNrNeg;
-
-		//! size of a single set of distribution functions 
-		STCON int    cDfNum;
-		//! direction vector contain vecs for all spatial dirs, even if not used for LBM model
-		STCON int    cDirNum;
-
-		//! distribution functions directions 
-		typedef enum {
-			 cDirInv=  -1,
-			 cDirC  =  0,
-			 cDirN  =  1,
-			 cDirS  =  2,
-			 cDirE  =  3,
-			 cDirW  =  4,
-			 cDirT  =  5,
-			 cDirB  =  6,
-			 cDirNE =  7,
-			 cDirNW =  8,
-			 cDirSE =  9,
-			 cDirSW = 10,
-			 cDirNT = 11,
-			 cDirNB = 12,
-			 cDirST = 13,
-			 cDirSB = 14,
-			 cDirET = 15,
-			 cDirEB = 16,
-			 cDirWT = 17,
-			 cDirWB = 18
-		} dfDir;
-
-		/* Vector Order 3D:
-		 *  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
-		 *  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
-		 *  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
-		 *  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1
-		 */
-
-		/*! name of the dist. function 
-			 only for nicer output */
-		STCON char* dfString[ 19 ];
-
-		/*! index of normal dist func, not used so far?... */
-		STCON int dfNorm[ 19 ];
-
-		/*! index of inverse dist func, not fast, but useful... */
-		STCON int dfInv[ 19 ];
-
-		/*! index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefX[ 19 ];
-		/*! index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefY[ 19 ];
-		/*! index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefZ[ 19 ];
-
-		/*! dist func vectors */
-		STCON int dfVecX[ 27 ];
-		STCON int dfVecY[ 27 ];
-		STCON int dfVecZ[ 27 ];
-
-		/*! arrays as before with doubles */
-		STCON LbmFloat dfDvecX[ 27 ];
-		STCON LbmFloat dfDvecY[ 27 ];
-		STCON LbmFloat dfDvecZ[ 27 ];
-
-		/*! principal directions */
-		STCON int princDirX[ 2*3 ];
-		STCON int princDirY[ 2*3 ];
-		STCON int princDirZ[ 2*3 ];
-
-		/*! vector lengths */
-		STCON LbmFloat dfLength[ 19 ];
-
-		/*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
-		static LbmFloat dfEquil[ 19 ];
-
-		/*! arrays for les model coefficients */
-		static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
-		static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
-
-}; // LbmData3D
-
-
-
-/*****************************************************************************/
-//! class for solver templating - 2D implementation 
-//class LbmD2Q9 : public LbmSolverInterface {
-class LbmD2Q9 {
-	
-	public:
-
-		// constructor, init interface
-		LbmD2Q9() {};
-		// virtual destructor 
-		virtual ~LbmD2Q9() {};
-		//! id string of solver
-		string getIdString() { return string("2D"); }
-
-		//! how many dimensions?
-		STCON int cDimension;
-
-		//! Wi factors for collide step 
-		STCON LbmFloat cCollenZero;
-		STCON LbmFloat cCollenOne;
-		STCON LbmFloat cCollenSqrtTwo;
-
-		//! threshold value for filled/emptied cells 
-		STCON LbmFloat cMagicNr2;
-		STCON LbmFloat cMagicNr2Neg;
-		STCON LbmFloat cMagicNr;
-		STCON LbmFloat cMagicNrNeg;
-
-		//! size of a single set of distribution functions 
-		STCON int    cDfNum;
-		STCON int    cDirNum;
-
-		//! distribution functions directions 
-		typedef enum {
-			 cDirInv=  -1,
-			 cDirC  =  0,
-			 cDirN  =  1,
-			 cDirS  =  2,
-			 cDirE  =  3,
-			 cDirW  =  4,
-			 cDirNE =  5,
-			 cDirNW =  6,
-			 cDirSE =  7,
-			 cDirSW =  8
-		} dfDir;
-
-		/* Vector Order 2D:
-		 * 0  1 2  3  4  5  6 7  8
-		 * 0, 0,0, 1,-1, 1,-1,1,-1 
-		 * 0, 1,-1, 0,0, 1,1,-1,-1  */
-
-		/* name of the dist. function 
-			 only for nicer output */
-		STCON char* dfString[ 9 ];
-
-		/* index of normal dist func, not used so far?... */
-		STCON int dfNorm[ 9 ];
-
-		/* index of inverse dist func, not fast, but useful... */
-		STCON int dfInv[ 9 ];
-
-		/* index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefX[ 9 ];
-		/* index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefY[ 9 ];
-		/* index of x reflected dist func for free slip, not valid for all DFs... */
-		STCON int dfRefZ[ 9 ];
-
-		/* dist func vectors */
-		STCON int dfVecX[ 9 ];
-		STCON int dfVecY[ 9 ];
-		/* Z, 2D values are all 0! */
-		STCON int dfVecZ[ 9 ];
-
-		/* arrays as before with doubles */
-		STCON LbmFloat dfDvecX[ 9 ];
-		STCON LbmFloat dfDvecY[ 9 ];
-		/* Z, 2D values are all 0! */
-		STCON LbmFloat dfDvecZ[ 9 ];
-
-		/*! principal directions */
-		STCON int princDirX[ 2*2 ];
-		STCON int princDirY[ 2*2 ];
-		STCON int princDirZ[ 2*2 ];
-
-		/* vector lengths */
-		STCON LbmFloat dfLength[ 9 ];
-
-		/* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
-		static LbmFloat dfEquil[ 9 ];
-
-		/*! arrays for les model coefficients */
-		static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
-		static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
-
-}; // LbmData3D
-
-
-
-// not needed hereafter
-#undef STCON
-
-
-
-/*****************************************************************************/
-//! class for solver templating - lbgk (srt) model implementation 
-template<class DQ>
-class LbmModelLBGK : public DQ , public LbmSolverInterface {
-	public:
-
-		/*! type for cells contents, needed for cell id interface */
-		typedef DQ LbmCellContents;
-		/*! type for cells */
-		typedef LbmCellTemplate< LbmCellContents > LbmCell;
-
-		// constructor
-		LbmModelLBGK() : DQ(), LbmSolverInterface() {};
-		// virtual destructor 
-		virtual ~LbmModelLBGK() {};
-		//! id string of solver
-		string getIdString() { return DQ::getIdString() + string("lbgk]"); }
-
-		/*! calculate length of velocity vector */
-		static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) {
-			return ((ux)*DQ::dfDvecX[l]+(uy)*DQ::dfDvecY[l]+(uz)*DQ::dfDvecZ[l]);
-		};
-
-		/*! calculate equilibrium DF for given values */
-		static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz) {
-			LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 
-			return( DQ::dfLength[l] *( 
-						+ rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
-						+ 3.0 *tmp 
-						+ 9.0/2.0 *(tmp*tmp) ) 
-					);
-		};
-
-
-		// input mux etc. as acceleration
-		// outputs rho,ux,uy,uz
-		/*inline void collideArrays_org(LbmFloat df[19], 				
-				LbmFloat &outrho, // out only!
-				// velocity modifiers (returns actual velocity!)
-				LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
-				LbmFloat omega
-			) {
-			LbmFloat rho=df[0]; 
-			LbmFloat ux = mux;
-			LbmFloat uy = muy;
-			LbmFloat uz = muz;
-			for(int l=1; l<DQ::cDfNum; l++) { 
-				rho += df[l]; 
-				ux  += (DQ::dfDvecX[l]*df[l]); 
-				uy  += (DQ::dfDvecY[l]*df[l]);  
-				uz  += (DQ::dfDvecZ[l]*df[l]);  
-			}  
-			for(int l=0; l<DQ::cDfNum; l++) { 
-				//LbmFloat tmp = (ux*DQ::dfDvecX[l]+uy*DQ::dfDvecY[l]+uz*DQ::dfDvecZ[l]); 
-				df[l] = (1.0-omega ) * df[l] + omega * ( getCollideEq(l,rho,ux,uy,uz) ); 
-			}  
-
-			mux = ux;
-			muy = uy;
-			muz = uz;
-			outrho = rho;
-		};*/
-		
-		// LES functions
-		inline LbmFloat getLesNoneqTensorCoeff(
-				LbmFloat df[], 				
-				LbmFloat feq[] ) {
-			LbmFloat Qo = 0.0;
-			for(int m=0; m< ((DQ::cDimension*DQ::cDimension)-DQ::cDimension)/2 ; m++) { 
-				LbmFloat qadd = 0.0;
-				for(int l=1; l<DQ::cDfNum; l++) { 
-					if(DQ::lesCoeffOffdiag[m][l]==0.0) continue;
-					qadd += DQ::lesCoeffOffdiag[m][l]*(df[l]-feq[l]);
-				}
-				Qo += (qadd*qadd);
-			}
-			Qo *= 2.0; // off diag twice
-			for(int m=0; m<DQ::cDimension; m++) { 
-				LbmFloat qadd = 0.0;
-				for(int l=1; l<DQ::cDfNum; l++) { 
-					if(DQ::lesCoeffDiag[m][l]==0.0) continue;
-					qadd += DQ::lesCoeffDiag[m][l]*(df[l]-feq[l]);
-				}
-				Qo += (qadd*qadd);
-			}
-			Qo = sqrt(Qo);
-			return Qo;
-		}
-		inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo) {
-			const LbmFloat tau = 1.0/omega;
-			const LbmFloat nu = (2.0*tau-1.0) * (1.0/6.0);
-			const LbmFloat C = csmago;
-			const LbmFloat Csqr = C*C;
-			LbmFloat S = -nu + sqrt( nu*nu + 18.0*Csqr*Qo ) / (6.0*Csqr);
-			return( 1.0/( 3.0*( nu+Csqr*S ) +0.5 ) );
-		}
-
-		// "normal" collision
-		inline void collideArrays(LbmFloat df[], 				
-				LbmFloat &outrho, // out only!
-				// velocity modifiers (returns actual velocity!)
-				LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
-				LbmFloat omega, LbmFloat csmago, LbmFloat *newOmegaRet = NULL
-			) {
-			LbmFloat rho=df[0]; 
-			LbmFloat ux = mux;
-			LbmFloat uy = muy;
-			LbmFloat uz = muz; 
-			for(int l=1; l<DQ::cDfNum; l++) { 
-				rho += df[l]; 
-				ux  += (DQ::dfDvecX[l]*df[l]); 
-				uy  += (DQ::dfDvecY[l]*df[l]);  
-				uz  += (DQ::dfDvecZ[l]*df[l]);  
-			}  
-			LbmFloat feq[19];
-			for(int l=0; l<DQ::cDfNum; l++) { 
-				feq[l] = getCollideEq(l,rho,ux,uy,uz); 
-			}
-
-			LbmFloat omegaNew;
-			if(csmago>0.0) {
-				LbmFloat Qo = getLesNoneqTensorCoeff(df,feq);
-				omegaNew = getLesOmega(omega,csmago,Qo);
-			} else {
-				omegaNew = omega; // smago off...
-			}
-			if(newOmegaRet) *newOmegaRet=omegaNew; // return value for stats
-
-			for(int l=0; l<DQ::cDfNum; l++) { 
-				df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
-			}  
-
-			mux = ux;
-			muy = uy;
-			muz = uz;
-			outrho = rho;
-		};
-
-}; // LBGK
-
-#ifdef LBMMODEL_DEFINED
-// force compiler error!
-ERROR - Dont include several LBM models at once...
-#endif
-#define LBMMODEL_DEFINED 1
-
-
-typedef LbmModelLBGK<  LbmD2Q9 > LbmBGK2D;
-typedef LbmModelLBGK< LbmD3Q19 > LbmBGK3D;
-
-
-#define LBMHEADER_H
-#endif
-
-
-

intern/elbeem/intern/lbmfunctions.h

@@ -1,320 +0,0 @@
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * All code distributed as part of El'Beem is covered by the version 2 of the 
- * GNU General Public License. See the file COPYING for details.
- * Copyright 2003-2005 Nils Thuerey
- *
- * Combined 2D/3D Lattice Boltzmann Solver templated helper functions
- * 
- *****************************************************************************/
-#ifndef LBMFUNCTIONS_H
-
-
-
-#if LBM_USE_GUI==1
-#define USE_GLUTILITIES
-#include "../gui/gui_utilities.h"
-
-//! display a single node
-template<typename D> 
-void 
-debugDisplayNode(fluidDispSettings *dispset, D *lbm, typename D::CellIdentifier cell ) {
-	//debugOut(" DD: "<<cell->getAsString() , 10);
-	ntlVec3Gfx org      = lbm->getCellOrigin( cell );
-	ntlVec3Gfx halfsize = lbm->getCellSize( cell );
-	int    set      = lbm->getCellSet( cell );
-	//debugOut(" DD: "<<cell->getAsString()<<" "<< (dispset->type) , 10);
-
-	bool     showcell = true;
-	int      linewidth = 1;
-	ntlColor col(0.5);
-	LbmFloat cscale = dispset->scale;
-
-#define DRAWDISPCUBE(col,scale) \
-	{	glLineWidth( linewidth ); \
-	  glColor3f( (col)[0], (col)[1], (col)[2]); \
-		ntlVec3Gfx s = org-(halfsize * (scale)); \
-		ntlVec3Gfx e = org+(halfsize * (scale)); \
-		drawCubeWire( s,e ); }
-
-	switch(dispset->type) {
-		case FLUIDDISPNothing: {
-				showcell = false;
-			} break;
-		case FLUIDDISPCelltypes: {
-				CellFlagType flag = lbm->getCellFlag(cell, set );
-				cscale = 0.5;
-
-				if(flag& CFInvalid  ) { if(!guiShowInvalid  ) return; }
-				if(flag& CFUnused   ) { if(!guiShowInvalid  ) return; }
-				if(flag& CFEmpty    ) { if(!guiShowEmpty    ) return; }
-				if(flag& CFInter    ) { if(!guiShowInterface) return; }
-				if(flag& CFNoDelete ) { if(!guiShowNoDelete ) return; }
-				if(flag& CFBnd      ) { if(!guiShowBnd      ) return; }
-
-				// only dismiss one of these types 
- 				if(flag& CFGrFromCoarse)  { if(!guiShowCoarseInner  ) return; } // inner not really interesting
-				else
-				if(flag& CFGrFromFine) { if(!guiShowCoarseBorder ) return; }
-				else
-				if(flag& CFFluid    )    { if(!guiShowFluid    ) return; }
-
-				if(flag& CFNoDelete) { // debug, mark nodel cells
-					ntlColor ccol(0.7,0.0,0.0);
-					DRAWDISPCUBE(ccol, 0.1);
-				}
-				if(flag& CFPersistMask) { // mark persistent flags
-					ntlColor ccol(0.5);
-					DRAWDISPCUBE(ccol, 0.125);
-				}
-				if(flag& CFNoBndFluid) { // mark persistent flags
-					ntlColor ccol(0,0,1);
-					DRAWDISPCUBE(ccol, 0.075);
-				}
-
-				/*if(flag& CFAccelerator) {
-					cscale = 0.55;
-					col = ntlColor(0,1,0);
-				} */
-				if(flag& CFInvalid) {
-					cscale = 0.50;
-					col = ntlColor(0.0,0,0.0);
-				}
-				/*else if(flag& CFSpeedSet) {
-					cscale = 0.55;
-					col = ntlColor(0.2,1,0.2);
-				}*/
-				else if(flag& CFBnd) {
-					cscale = 0.59;
-					col = ntlColor(0.4);
-				}
-
-				else if(flag& CFInter) {
-					cscale = 0.55;
-					col = ntlColor(0,1,1);
-
-				} else if(flag& CFGrFromCoarse) {
-					// draw as - with marker
-					ntlColor col2(0.0,1.0,0.3);
-					DRAWDISPCUBE(col2, 0.1);
-					cscale = 0.5;
-					showcell=false; // DEBUG
-				}
-				else if(flag& CFFluid) {
-					cscale = 0.5;
-					if(flag& CFGrToFine) {
-						ntlColor col2(0.5,0.0,0.5);
-						DRAWDISPCUBE(col2, 0.1);
-						col = ntlColor(0,0,1);
-					}
-					if(flag& CFGrFromFine) {
-						ntlColor col2(1.0,1.0,0.0);
-						DRAWDISPCUBE(col2, 0.1);
-						col = ntlColor(0,0,1);
-					} else if(flag& CFGrFromCoarse) {
-						// draw as fluid with marker
-						ntlColor col2(0.0,1.0,0.3);
-						DRAWDISPCUBE(col2, 0.1);
-						col = ntlColor(0,0,1);
-					} else {
-						col = ntlColor(0,0,1);
-					}
-				}
-				else if(flag& CFEmpty) {
-					showcell=false;
-				}
-
-			} break;
-		case FLUIDDISPVelocities: {
-				// dont use cube display
-				LbmVec vel = lbm->getCellVelocity( cell, set );
-				glBegin(GL_LINES);
-				glColor3f( 0.0,0.0,0.0 );
-				glVertex3f( org[0], org[1], org[2] );
-				org += vec2G(vel * 10.0 * cscale);
-				glColor3f( 1.0,1.0,1.0 );
-				glVertex3f( org[0], org[1], org[2] );
-				glEnd();
-				showcell = false;
-			} break;
-		case FLUIDDISPCellfills: {
-				CellFlagType flag = lbm->getCellFlag( cell,set );
-				cscale = 0.5;
-
-				if(flag& CFFluid) {
-					cscale = 0.75;
-					col = ntlColor(0,0,0.5);
-				}
-				else if(flag& CFInter) {
-					cscale = 0.75 * lbm->getCellMass(cell,set);
-					col = ntlColor(0,1,1);
-				}
-				else {
-					showcell=false;
-				}
-
-					if( ABS(lbm->getCellMass(cell,set)) < 10.0 ) {
-						cscale = 0.75 * lbm->getCellMass(cell,set);
-					} else {
-						showcell = false;
-					}
-					if(cscale>0.0) {
-						col = ntlColor(0,1,1);
-					} else {
-						col = ntlColor(1,1,0);
-					}
-			// TODO
-			} break;
-		case FLUIDDISPDensity: {
-				LbmFloat rho = lbm->getCellDensity(cell,set);
-				cscale = rho*rho * 0.25;
-				col = ntlColor( MIN(0.5+cscale,1.0) , MIN(0.0+cscale,1.0), MIN(0.0+cscale,1.0) );
-				cscale *= 2.0;
-			} break;
-		case FLUIDDISPGrid: {
-				cscale = 0.59;
-				col = ntlColor(1.0);
-			} break;
-		default: {
-				cscale = 0.5;
-				col = ntlColor(1.0,0.0,0.0);
-			} break;
-	}
-
-	if(!showcell) return;
-	DRAWDISPCUBE(col, cscale);
-}
-
-//! debug display function
-//  D has to implement the CellIterator interface
-template<typename D>
-void 
-lbmDebugDisplay(fluidDispSettings *dispset, D *lbm) {
-	//je nach solver...?
-	if(!dispset->on) return;
-	glDisable( GL_LIGHTING ); // dont light lines
-
-	typename D::CellIdentifier cid = lbm->getFirstCell();
-	for(; lbm->noEndCell( cid );
-	      lbm->advanceCell( cid ) ) {
-		// display...
-#if (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)
-		::debugDisplayNode<>(dispset, lbm, cid );
-#else
-		debugDisplayNode<D>(dispset, lbm, cid );
-#endif
-	}
-	delete cid;
-
-	glEnable( GL_LIGHTING ); // dont light lines
-}
-
-//! debug display function
-//  D has to implement the CellIterator interface
-template<typename D>
-void 
-lbmMarkedCellDisplay(D *lbm) {
-	fluidDispSettings dispset;
-	// trick - display marked cells as grid displa -> white, big
-	dispset.type = FLUIDDISPGrid;
-	dispset.on = true;
-	glDisable( GL_LIGHTING ); // dont light lines
-	
-	typename D::CellIdentifier cid = lbm->markedGetFirstCell();
-	while(cid) {
-	//for(; lbm->markedNoEndCell( cid );
-	      //cid = lbm->markedAdvanceCell( cid ) ) {
-		// display... FIXME? this is a bit inconvenient...
-		//MarkedCellIdentifier *mid = dynamic_cast<MarkedCellIdentifier *>( cid );
-#if (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)
-		//::debugDisplayNode<>(&dispset, lbm, mid->mpCell );
-		::debugDisplayNode<>(&dispset, lbm, cid );
-#else
-		//debugDisplayNode<D>(&dispset, lbm, mid->mpCell );
-		debugDisplayNode<D>(&dispset, lbm, cid );
-#endif
-		cid = lbm->markedAdvanceCell();
-	}
-	delete cid;
-
-	glEnable( GL_LIGHTING ); // dont light lines
-}
-
-#endif // LBM_USE_GUI
-
-//! display a single node
-template<typename D> 
-void 
-debugPrintNodeInfo(D *lbm, typename D::CellIdentifier cell, string printInfo,
-		// force printing of one set? default = -1 = off
-		int forceSet=-1) {
-  bool printDF     = false;
-  bool printRho    = false;
-  bool printVel    = false;
-  bool printFlag   = false;
-  bool printGeom   = false;
-  bool printMass=false;
-	bool printBothSets = false;
-
-	for(size_t i=0; i<printInfo.length()-0; i++) {
-		char what = printInfo[i];
-		switch(what) {
-			case '+': // all on
-								printDF = true; printRho = true; printVel = true; printFlag = true; printGeom = true; printMass = true; 
-								printBothSets = true; break;
-			case '-': // all off
-								printDF = false; printRho = false; printVel = false; printFlag = false; printGeom = false; printMass = false; 
-								printBothSets = false; break;
-			case 'd': printDF = true; break;
-			case 'r': printRho = true; break;
-			case 'v': printVel = true; break;
-			case 'f': printFlag = true; break;
-			case 'g': printGeom = true; break;
-			case 'm': printMass = true; break;
-			case 's': printBothSets = true; break;
-			default: 
-				errFatal("debugPrintNodeInfo","Invalid node info id "<<what,SIMWORLD_GENERICERROR); return;
-		}
-	}
-
-	ntlVec3Gfx org      = lbm->getCellOrigin( cell );
-	ntlVec3Gfx halfsize = lbm->getCellSize( cell );
-	int    set      = lbm->getCellSet( cell );
-	debMsgStd("debugPrintNodeInfo",DM_NOTIFY, "Printing cell info '"<<printInfo<<"' for node: "<<cell->getAsString()<<" from "<<lbm->getName()<<" currSet:"<<set , 1);
-	if(printGeom) debMsgStd("                  ",DM_MSG, "Org:"<<org<<" Halfsize:"<<halfsize<<" ", 1);
-
-	int setmax = 2;
-	if(!printBothSets) setmax = 1;
-	if(forceSet>=0) setmax = 1;
-
-	for(int s=0; s<setmax; s++) {
-		int workset = set;
-		if(s==1){ workset = (set^1); }		
-		if(forceSet>=0) workset = forceSet;
-		debMsgStd("                  ",DM_MSG, "Printing set:"<<workset<<" orgSet:"<<set, 1);
-		
-		if(printDF) {
-			for(int l=0; l<lbm->getDfNum(); l++) { // FIXME ??
-				debMsgStd("                  ",DM_MSG, "  Df"<<l<<": "<<lbm->getCellDf(cell,workset,l), 1);
-			}
-		}
-		if(printRho) {
-			debMsgStd("                  ",DM_MSG, "  Rho: "<<lbm->getCellDensity(cell,workset), 1);
-		}
-		if(printVel) {
-			debMsgStd("                  ",DM_MSG, "  Vel: "<<lbm->getCellVelocity(cell,workset), 1);
-		}
-		if(printFlag) {
-			CellFlagType flag = lbm->getCellFlag(cell,workset);
-			debMsgStd("                  ",DM_MSG, "  Flg: "<< flag<<" "<<convertFlags2String( flag ) <<" "<<convertCellFlagType2String( flag ), 1);
-		}
-		if(printMass) {
-			debMsgStd("                  ",DM_MSG, "  Mss: "<<lbm->getCellMass(cell,workset), 1);
-		}
-	}
-}
-
-#define LBMFUNCTIONS_H
-#endif
-

intern/elbeem/intern/ntl_blenderdumper.cpp

@@ -24,7 +24,7 @@
  * Constructor
  *****************************************************************************/
 ntlBlenderDumper::ntlBlenderDumper(string filename, bool commandlineMode) :
-	ntlRaytracer(filename,commandlineMode),
+	ntlWorld(filename,commandlineMode),
 	mpTrafo(NULL)
 {
   ntlRenderGlobals *glob = mpGlob;
@@ -194,7 +194,7 @@ int ntlBlenderDumper::renderScene( void )
 	// still render for preview...
 	if(debugOut) {
 		debMsgStd("ntlBlenderDumper::renderScene",DM_NOTIFY,"Performing preliminary render", 1);
-		ntlRaytracer::renderScene(); }
+		ntlWorld::renderScene(); }
 	else {
 		// next frame 
 		glob->setAniCount( glob->getAniCount() +1 );

intern/elbeem/intern/ntl_blenderdumper.h

@@ -7,12 +7,12 @@
  *
  *****************************************************************************/
 #ifndef NTL_BLENDERDUMPER_H
-#include "ntl_raytracer.h"
+#include "ntl_world.h"
 
 template<class Scalar> class ntlMatrix4x4;
 
 class ntlBlenderDumper :
-	public ntlRaytracer
+	public ntlWorld
 {
 public:
   /*! Constructor */

intern/elbeem/intern/ntl_bsptree.cpp

@@ -579,18 +579,6 @@ void ntlTree::intersect(const ntlRay &ray, gfxReal &distance,
 						mint = t;	  
 						hit = (*iter);
 						mintu = u; mintv = v;
-						
-						if((ray.getRenderglobals())&&(ray.getRenderglobals()->getDebugOut() > 5)) {  // DEBUG!!!
-							errorOut("Tree tri hit at "<<t<<","<<mint<<" triangle: "<<PRINT_TRIANGLE( (*hit), (*mpVertices) ) );
-							gfxReal u1=0.0,v1=0.0, t1=-1.0;
-							ray.intersectTriangle( mpVertices, hit, t1,u1,v1);
-							errorOut("Tree second test1 :"<<t1<<" u1:"<<u1<<" v1:"<<v1 );
-							if(t==GFX_REAL_MAX) errorOut( "Tree MAX t " );
-							//errorOut( mpVertices[ (*iter).getPoints()[0] ][0] );
-						}
-
-						//retnormal = -(e2-e0).crossProd(e1-e0); // DEBUG
-
 					}
 				}
 

intern/elbeem/intern/ntl_geometryobject.cpp

@@ -85,6 +85,7 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
 		mGeoInitId = -1;
 	}
 	mInitialVelocity = vec2G( mpAttrs->readVec3d("initial_velocity", vec2D(mInitialVelocity),"ntlGeometryObject", "mInitialVelocity", false));
+	mGeoPartSlipValue = mpAttrs->readFloat("geoinit_partslip", mGeoPartSlipValue,"ntlGeometryObject", "mGeoPartSlipValue", false);
 	debMsgStd("ntlGeometryObject::initialize",DM_MSG,"GeoObj '"<<this->getName()<<"': gid="<<mGeoInitId<<" gtype="<<mGeoInitType<<","<<ginitStr<<
 			" gvel="<<mInitialVelocity<<" gisect="<<mGeoInitIntersect, 10); // debug
 

intern/elbeem/intern/ntl_geometryobject.h

@@ -72,6 +72,10 @@ class ntlGeometryObject : public ntlGeometryClass
 		inline bool getGeoInitIntersect() const { return mGeoInitIntersect; }
 		inline void setGeoInitIntersect(bool set) { mGeoInitIntersect=set; }
 
+		/*! Set/get the part slip value*/
+		inline bool getGeoPartSlipValue() const { return mGeoPartSlipValue; }
+		inline void setGeoPartSlipValue(float set) { mGeoPartSlipValue=set; }
+
 	protected:
 
 		/*! Point to a property object describing the surface of this object */
@@ -94,6 +98,8 @@ class ntlGeometryObject : public ntlGeometryClass
 		ntlVec3Gfx mInitialVelocity;
 		/*! perform more accurate intersecting geo init for this object? */
 		bool mGeoInitIntersect;
+		/*! part slip bc value */
+		float mGeoPartSlipValue;
 
 	public:
 

intern/elbeem/intern/ntl_ray.cpp

@@ -466,7 +466,7 @@ const ntlColor ntlRay::shade() //const
   if (closest != NULL) {
 
 		ntlVec3Gfx triangleNormal = tri->getNormal();
-		if( equal(triangleNormal, ntlVec3Gfx(0.0)) ) errorOut("ntlRaytracer warning: trinagle normal= 0 "); // DEBUG
+		if( equal(triangleNormal, ntlVec3Gfx(0.0)) ) errorOut("ntlRay warning: trinagle normal= 0 "); // DEBUG
 		/* intersection on inside faces? if yes invert normal afterwards */
 		gfxReal valDN; // = mDirection | normal;
 		valDN = dot(mDirection, triangleNormal);

intern/elbeem/intern/ntl_world.cpp

@@ -11,7 +11,7 @@
 #include <sys/stat.h>
 #include <sstream>
 #include "utilities.h"
-#include "ntl_raytracer.h"
+#include "ntl_world.h"
 #include "ntl_scene.h"
 #include "parametrizer.h"
 #include "globals.h"
@@ -34,7 +34,7 @@ void setPointers( ntlRenderGlobals *setglob);
 /******************************************************************************
  * Constructor
  *****************************************************************************/
-ntlRaytracer::ntlRaytracer(string filename, bool commandlineMode) :
+ntlWorld::ntlWorld(string filename, bool commandlineMode) :
 	mpGlob(NULL), 
   mpLightList(NULL), mpPropList(NULL), mpSims(NULL),
 	mpOpenGLRenderer(NULL),
@@ -71,7 +71,7 @@ ntlRaytracer::ntlRaytracer(string filename, bool commandlineMode) :
   long sstartTime = getTime();
 	scene->buildScene();
 	long sstopTime = getTime();
-	debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"Scene build time: "<< getTimeString(sstopTime-sstartTime) <<" ", 10);
+	debMsgStd("ntlWorld::ntlWorld",DM_MSG,"Scene build time: "<< getTimeString(sstopTime-sstartTime) <<" ", 10);
 
 	// TODO check simulations, run first steps
 	mFirstSim = -1;
@@ -89,31 +89,31 @@ ntlRaytracer::ntlRaytracer(string filename, bool commandlineMode) :
 			if(mFirstSim>=0) {
 				if( (*mpSims)[i]->getStepTime() > (*mpSims)[mFirstSim]->getStepTime() ) {
 					mFirstSim = i;
-					debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"First Sim changed: "<<i ,10);
+					debMsgStd("ntlWorld::ntlWorld",DM_MSG,"First Sim changed: "<<i ,10);
 				}
 			}
 			// check any valid sim
 			if(mFirstSim<0) {
 				mFirstSim = i;
-				debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"First Sim: "<<i ,10);
+				debMsgStd("ntlWorld::ntlWorld",DM_MSG,"First Sim: "<<i ,10);
 			}
 		}
 
 		if(mFirstSim>=0) {
-			debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"Anistart Time: "<<(*mpSims)[mFirstSim]->getStartTime() ,10);
+			debMsgStd("ntlWorld::ntlWorld",DM_MSG,"Anistart Time: "<<(*mpSims)[mFirstSim]->getStartTime() ,10);
 			while(mSimulationTime < (*mpSims)[mFirstSim]->getStartTime() ) {
-			debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"Anistart Time: "<<(*mpSims)[mFirstSim]->getStartTime()<<" simtime:"<<mSimulationTime ,10);
+			debMsgStd("ntlWorld::ntlWorld",DM_MSG,"Anistart Time: "<<(*mpSims)[mFirstSim]->getStartTime()<<" simtime:"<<mSimulationTime ,10);
 				advanceSims();
 			}
 			long stopTime = getTime();
 
 			mSimulationTime += (*mpSims)[mFirstSim]->getStartTime();
-			debMsgStd("ntlRaytracer::ntlRaytracer",DM_MSG,"Time for start simulations times "<<": "<< getTimeString(stopTime-startTime) <<"s ", 1);
+			debMsgStd("ntlWorld::ntlWorld",DM_MSG,"Time for start simulations times "<<": "<< getTimeString(stopTime-startTime) <<"s ", 1);
 #ifndef NOGUI
 			guiResetSimulationTimeRange( mSimulationTime );
 #endif
 		} else {
-			if(!mpGlob->getSingleFrameMode()) debMsgStd("ntlRaytracer::ntlRaytracer",DM_WARNING,"No active simulations!", 1);
+			if(!mpGlob->getSingleFrameMode()) debMsgStd("ntlWorld::ntlWorld",DM_WARNING,"No active simulations!", 1);
 		}
 	}
 
@@ -132,7 +132,7 @@ ntlRaytracer::ntlRaytracer(string filename, bool commandlineMode) :
 /******************************************************************************
  * Destructor
  *****************************************************************************/
-ntlRaytracer::~ntlRaytracer()
+ntlWorld::~ntlWorld()
 {
 	delete mpGlob->getScene();
   delete mpGlob;
@@ -146,7 +146,7 @@ ntlRaytracer::~ntlRaytracer()
 
 /******************************************************************************/
 /*! set single frame rendering to filename */
-void ntlRaytracer::setSingleFrameOut(string singleframeFilename) {
+void ntlWorld::setSingleFrameOut(string singleframeFilename) {
 	mpGlob->setSingleFrameMode(true);
 	mpGlob->setSingleFrameFilename(singleframeFilename);
 }
@@ -159,17 +159,17 @@ extern "C" {
 	void simulateThreadIncreaseFrame(void);
 }
 #endif // ELBEEM_BLENDER==1
-int ntlRaytracer::renderAnimation( void )
+int ntlWorld::renderAnimation( void )
 {
 	// only single pic currently
-	//debMsgStd("ntlRaytracer::renderAnimation : Warning only simulating...",1);
+	//debMsgStd("ntlWorld::renderAnimation : Warning only simulating...",1);
  	if(mpGlob->getAniFrames() < 0) {
-		debMsgStd("ntlRaytracer::renderAnimation",DM_NOTIFY,"No frames to render... ",1);
+		debMsgStd("ntlWorld::renderAnimation",DM_NOTIFY,"No frames to render... ",1);
 		return 1;
 	}
 
 	if(mFirstSim<0) {
-		debMsgStd("ntlRaytracer::renderAnimation",DM_NOTIFY,"No reference animation found...",1);
+		debMsgStd("ntlWorld::renderAnimation",DM_NOTIFY,"No reference animation found...",1);
 		return 1;
 	} 
 
@@ -184,7 +184,7 @@ int ntlRaytracer::renderAnimation( void )
 #endif // ELBEEM_BLENDER==1
 
 		if(mpSims->size() <= 0) {
-			debMsgStd("ntlRaytracer::renderAnimation",DM_NOTIFY,"No simulations found, stopping...",1);
+			debMsgStd("ntlWorld::renderAnimation",DM_NOTIFY,"No simulations found, stopping...",1);
 			return 1;
 		}
 	}
@@ -198,7 +198,7 @@ int ntlRaytracer::renderAnimation( void )
  * this function is run in another thread, and communicates 
  * with the parent thread via a mutex 
  *****************************************************************************/
-int ntlRaytracer::renderVisualization( bool multiThreaded ) 
+int ntlWorld::renderVisualization( bool multiThreaded ) 
 {
 #ifndef NOGUI
 	//gfxReal deltat = 0.0015;
@@ -206,7 +206,7 @@ int ntlRaytracer::renderVisualization( bool multiThreaded )
 	while(!getStopRenderVisualization()) {
 
 		if(mpSims->size() <= 0) {
-			debMsgStd("ntlRaytracer::renderVisualization",DM_NOTIFY,"No simulations found, stopping...",1);
+			debMsgStd("ntlWorld::renderVisualization",DM_NOTIFY,"No simulations found, stopping...",1);
 			stopSimulationThread();
 			break;
 		}
@@ -216,7 +216,7 @@ int ntlRaytracer::renderVisualization( bool multiThreaded )
 			long startTime = getTime();
 			advanceSims();
 			long stopTime = getTime();
-			debMsgStd("ntlRaytracer::renderVisualization",DM_MSG,"Time for "<<mSimulationTime<<": "<< getTimeString(stopTime-startTime) <<"s ", 10);
+			debMsgStd("ntlWorld::renderVisualization",DM_MSG,"Time for "<<mSimulationTime<<": "<< getTimeString(stopTime-startTime) <<"s ", 10);
 		} else {
 			double targetTime = mSimulationTime + (*mpSims)[mFirstSim]->getStepTime();
 			singleStepSims(targetTime);
@@ -227,11 +227,11 @@ int ntlRaytracer::renderVisualization( bool multiThreaded )
 				if(!(*mpSims)[i]->getPanic()) allPanic = false;
 			}
 			if(allPanic) {
-				warnMsg("ntlRaytracer::advanceSims","All sims panicked... stopping thread" );
+				warnMsg("ntlWorld::advanceSims","All sims panicked... stopping thread" );
 				setStopRenderVisualization( true );
 			}
 			if(!SIMWORLD_OK()) {
-				warnMsg("ntlRaytracer::advanceSims","World state error... stopping" );
+				warnMsg("ntlWorld::advanceSims","World state error... stopping" );
 				setStopRenderVisualization( true );
 			}
 		}
@@ -254,7 +254,7 @@ int ntlRaytracer::renderVisualization( bool multiThreaded )
 	return 0;
 }
 /*! render a single step for viz mode */
-int ntlRaytracer::singleStepVisualization( void ) 
+int ntlWorld::singleStepVisualization( void ) 
 {
 	mThreadRunning = true;
 	double targetTime = mSimulationTime + (*mpSims)[mFirstSim]->getStepTime();
@@ -279,12 +279,12 @@ int ntlRaytracer::singleStepVisualization( void )
 /******************************************************************************
  * advance simulations by time t 
  *****************************************************************************/
-int ntlRaytracer::advanceSims()
+int ntlWorld::advanceSims()
 {
 	bool done = false;
 	bool allPanic = true;
 	double targetTime = mSimulationTime + (*mpSims)[mFirstSim]->getFrameTime();
-	//debMsgStd("ntlRaytracer::advanceSims",DM_MSG,"Advancing sims to "<<targetTime, 10 ); // timedebug
+	//debMsgStd("ntlWorld::advanceSims",DM_MSG,"Advancing sims to "<<targetTime, 10 ); // timedebug
 
 	while(!done) {
 		double nextTargetTime = (*mpSims)[mFirstSim]->getCurrentTime() + (*mpSims)[mFirstSim]->getStepTime();
@@ -296,14 +296,14 @@ int ntlRaytracer::advanceSims()
 		for(size_t i=0;i<mpSims->size();i++) {
 			if(!(*mpSims)[i]->getVisible()) continue;
 			if((*mpSims)[i]->getPanic()) allPanic = true; // do any panic now!?
-			//debMsgStd("ntlRaytracer::advanceSims",DM_MSG, " sim "<<i<<" c"<<(*mpSims)[i]->getCurrentTime()<<" p"<<(*mpSims)[i]->getPanic()<<" t"<<targetTime, 10); // debug // timedebug
+			//debMsgStd("ntlWorld::advanceSims",DM_MSG, " sim "<<i<<" c"<<(*mpSims)[i]->getCurrentTime()<<" p"<<(*mpSims)[i]->getPanic()<<" t"<<targetTime, 10); // debug // timedebug
 		}
 		if( (targetTime - (*mpSims)[mFirstSim]->getCurrentTime()) > LBM_TIME_EPSILON) done=false;
 		if(allPanic) done = true;
 	}
 
 	if(allPanic) {
-		warnMsg("ntlRaytracer::advanceSims","All sims panicked... stopping thread" );
+		warnMsg("ntlWorld::advanceSims","All sims panicked... stopping thread" );
 		setStopRenderVisualization( true );
 	}
 	for(size_t i=0;i<mpSims->size();i++) {
@@ -318,10 +318,10 @@ int ntlRaytracer::advanceSims()
 
 /* advance simulations by a single step */
 /* dont check target time, if *targetTime==NULL */
-void ntlRaytracer::singleStepSims(double targetTime) {
+void ntlWorld::singleStepSims(double targetTime) {
 	const bool debugTime = false;
 	//double targetTime = mSimulationTime + (*mpSims)[mFirstSim]->getStepTime();
-	if(debugTime) errMsg("ntlRaytracer::singleStepSims","Target time: "<<targetTime);
+	if(debugTime) errMsg("ntlWorld::singleStepSims","Target time: "<<targetTime);
 
 	for(size_t i=0;i<mpSims->size();i++) {
 		SimulationObject *sim = (*mpSims)[i];
@@ -330,9 +330,9 @@ void ntlRaytracer::singleStepSims(double targetTime) {
 		bool done = false;
 		while(!done) {
 			// try to prevent round off errs
-			if(debugTime) errMsg("ntlRaytracer::singleStepSims","Test sim "<<i<<" curt:"<< sim->getCurrentTime()<<" target:"<<targetTime<<" delta:"<<(targetTime - sim->getCurrentTime())<<" stept:"<<sim->getStepTime()<<" leps:"<<LBM_TIME_EPSILON ); // timedebug
+			if(debugTime) errMsg("ntlWorld::singleStepSims","Test sim "<<i<<" curt:"<< sim->getCurrentTime()<<" target:"<<targetTime<<" delta:"<<(targetTime - sim->getCurrentTime())<<" stept:"<<sim->getStepTime()<<" leps:"<<LBM_TIME_EPSILON ); // timedebug
 			if( (targetTime - sim->getCurrentTime()) > LBM_TIME_EPSILON) {
-				if(debugTime) errMsg("ntlRaytracer::singleStepSims","Stepping sim "<<i<<" t:"<< sim->getCurrentTime()); // timedebug
+				if(debugTime) errMsg("ntlWorld::singleStepSims","Stepping sim "<<i<<" t:"<< sim->getCurrentTime()); // timedebug
 				sim->step();
 			} else {
 				done = true;
@@ -353,7 +353,7 @@ void ntlRaytracer::singleStepSims(double targetTime) {
  * Render the current scene
  * uses the global variables from the parser
  *****************************************************************************/
-int ntlRaytracer::renderScene( void )
+int ntlWorld::renderScene( void )
 {
 #ifndef ELBEEM_BLENDER
 	char nrStr[5];														/* nr conversion */
@@ -374,7 +374,7 @@ int ntlRaytracer::renderScene( void )
 		if(mpGlob->getFrameSkip()) {
 			struct stat statBuf;
 			if(stat(outfn_conv.str().c_str(),&statBuf) == 0) {
-				errorOut("ntlRaytracer::renderscene Warning: file "<<outfn_conv.str()<<" already exists - skipping frame..."); 
+				errorOut("ntlWorld::renderscene Warning: file "<<outfn_conv.str()<<" already exists - skipping frame..."); 
 				glob->setAniCount( glob->getAniCount() +1 );
 				return(2);
 			}
@@ -472,7 +472,7 @@ int ntlRaytracer::renderScene( void )
 	glob->setCounterSceneInter(0);
 	for (int scanline=Yres ; scanline > 0 ; --scanline) {
     
-		debugOutInter( "ntlRaytracer::renderScene: Line "<<scanline<<
+		debugOutInter( "ntlWorld::renderScene: Line "<<scanline<<
 								 " ("<< ((Yres-scanline)*100/Yres) <<"%) ", 2, 2000 );
 		screenPos = glob->getLookat() + upVec*((2.0*scanline-Yres)/Yres)
 			- rightVec;
@@ -679,7 +679,7 @@ int ntlRaytracer::renderScene( void )
 #ifndef NOPNG
 		writePng(outfn_conv.str().c_str(), rows, w, h);
 #else // NOPNG
-		debMsgStd("ntlRaytracer::renderScene",DM_NOTIFY, "No PNG linked, no picture...", 1);
+		debMsgStd("ntlWorld::renderScene",DM_NOTIFY, "No PNG linked, no picture...", 1);
 #endif // NOPNG
 	}
 
@@ -696,7 +696,7 @@ int ntlRaytracer::renderScene( void )
 				 getTimeString(totalStart-timeStart).c_str(), getTimeString(rendEnd-rendStart).c_str(), getTimeString(timeEnd-timeStart).c_str(),
 				 glob->getCounterShades(),
 				 glob->getCounterSceneInter() );
-	debMsgStd("ntlRaytracer::renderScene",DM_MSG, resout, 1 );
+	debMsgStd("ntlWorld::renderScene",DM_MSG, resout, 1 );
 
 	/* clean stuff up */
 	delete [] aaCol;
@@ -705,7 +705,7 @@ int ntlRaytracer::renderScene( void )
 	glob->getScene()->cleanupScene();
 
 	if(mpGlob->getSingleFrameMode() ) {
-		debMsgStd("ntlRaytracer::renderScene",DM_NOTIFY, "Single frame mode done...", 1 );
+		debMsgStd("ntlWorld::renderScene",DM_NOTIFY, "Single frame mode done...", 1 );
 		return 1;
 	}
 #endif // ELBEEM_BLENDER

intern/elbeem/intern/ntl_world.h

@@ -18,13 +18,13 @@
 #include "simulation_object.h"
 class ntlOpenGLRenderer;
 
-class ntlRaytracer
+class ntlWorld
 {
 	public:
 		/*! Constructor */
-		ntlRaytracer(string filename, bool commandlineMode);
+		ntlWorld(string filename, bool commandlineMode);
 		/*! Destructor */
-		virtual ~ntlRaytracer( void );
+		virtual ~ntlWorld( void );
 
 		/*! render a whole animation (command line mode) */
 		int renderAnimation( void );

intern/elbeem/intern/simulation_object.cpp

@@ -10,8 +10,6 @@
 #include "simulation_object.h"
 #include "ntl_bsptree.h"
 #include "ntl_scene.h"
-#include "factory_lbm.h"
-#include "lbmfunctions.h"
 
 #ifdef _WIN32
 #else
@@ -19,6 +17,8 @@
 #endif
 
 
+//! lbm factory functions
+LbmSolverInterface* createSolver();
 
 
 /******************************************************************************
@@ -109,13 +109,8 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 	//mDimension is deprecated
 	mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false ); 
 	if(mSolverType == stnFsgr) {
-		mpLbm = createSolverLbmFsgr(); 
+		mpLbm = createSolver(); 
 	} else if(mSolverType == stnOld) {
-#ifdef LBM_INCLUDE_TESTSOLVERS
-		// old solver for gfx demo
-		mpLbm = createSolverOld();
-#endif // LBM_TESTSOLVER
-	} else {
 		errFatal("SimulationObject::initializeLbmSimulation","Invalid solver type - note that mDimension is deprecated, use the 'solver' keyword instead", SIMWORLD_INITERROR);
 		return 1;
 	}
@@ -311,7 +306,8 @@ void SimulationObject::drawDebugDisplay() {
 	//errorOut( mDebugType <<"//"<< mDebDispSet[mDebugType].type );
 	mpLbm->debugDisplay( &mDebDispSet[ mDebugType ] );
 
-	::lbmMarkedCellDisplay<>( mpLbm );
+	//::lbmMarkedCellDisplay<>( mpLbm );
+	mpLbm->lbmMarkedCellDisplay();
 #endif
 }
 
@@ -322,7 +318,7 @@ void SimulationObject::drawInteractiveDisplay()
 	if(!getVisible()) return;
 	if(mSelectedCid) {
 		// in debugDisplayNode if dispset is on is ignored...
-		::debugDisplayNode<>( &mDebDispSet[ FLUIDDISPGrid ], mpLbm, mSelectedCid );
+		mpLbm->debugDisplayNode( &mDebDispSet[ FLUIDDISPGrid ], mSelectedCid );
 	}
 #endif
 }
@@ -351,7 +347,8 @@ void SimulationObject::setMousePos(int x,int y, ntlVec3Gfx org, ntlVec3Gfx dir)
 void SimulationObject::setMouseClick()
 {
 	if(mSelectedCid) {
-		::debugPrintNodeInfo<>( mpLbm, mSelectedCid, mpLbm->getNodeInfoString() );
+		//::debugPrintNodeInfo<>( mpLbm, mSelectedCid, mpLbm->getNodeInfoString() );
+		mpLbm->debugPrintNodeInfo( mSelectedCid );
 	}
 }
 

intern/elbeem/intern/simulation_object.h

@@ -7,13 +7,13 @@
  *
  *****************************************************************************/
 
-#ifndef ELBEEM_SIMINTERFACE
+#ifndef SIMULATION_OBJECT_H
-#define ELBEEM_SIMINTERFACE
+#define SIMULATION_OBJECT_H
 
 
 #define USE_GLUTILITIES
 #include "ntl_geometryshader.h"
-#include "lbmdimensions.h"
+#include "solver_interface.h"
 #include "parametrizer.h"
 #include "particletracer.h"
 

intern/elbeem/intern/solver_class.h

@@ -0,0 +1,690 @@
+/******************************************************************************
+ *
+ * El'Beem - the visual lattice boltzmann freesurface simulator
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2005 Nils Thuerey
+ *
+ * Combined 2D/3D Lattice Boltzmann standard solver classes
+ *
+ *****************************************************************************/
+
+
+#ifndef LBMFSGRSOLVER_H
+
+// blender interface
+#if ELBEEM_BLENDER==1
+// warning - for MSVC this has to be included
+// _before_ ntl_vector3dim
+#include "SDL.h"
+#include "SDL_thread.h"
+#include "SDL_mutex.h"
+extern "C" {
+	void simulateThreadIncreaseFrame(void);
+	extern SDL_mutex *globalBakeLock;
+	extern int        globalBakeState;
+	extern int        globalBakeFrame;
+}
+#endif // ELBEEM_BLENDER==1
+
+#include "utilities.h"
+#include "solver_interface.h"
+#include "ntl_scene.h"
+#include <stdio.h>
+
+#if PARALLEL==1
+#include <omp.h>
+#endif // PARALLEL=1
+#ifndef PARALLEL
+#define PARALLEL 0
+#endif // PARALLEL
+
+#ifndef LBMMODEL_DEFINED
+// force compiler error!
+ERROR - define model first!
+#endif // LBMMODEL_DEFINED
+
+
+// general solver setting defines
+
+//! debug coordinate accesses and the like? (much slower)
+#define FSGR_STRICT_DEBUG 0
+
+//! debug coordinate accesses and the like? (much slower)
+#define FSGR_OMEGA_DEBUG 0
+
+//! OPT3D quick LES on/off, only debug/benchmarking
+#define USE_LES 1
+
+//! order of interpolation (0=always current/1=interpolate/2=always other)
+//#define TIMEINTORDER 0
+// TODO remove interpol t param, also interTime
+
+//! refinement border method (1 = small border / 2 = larger)
+#define REFINEMENTBORDER 1
+
+// use optimized 3D code?
+#if LBMDIM==2
+#define OPT3D 0
+#else
+// determine with debugging...
+#	if FSGR_STRICT_DEBUG==1
+#		define OPT3D 0
+#	else // FSGR_STRICT_DEBUG==1
+// usually switch optimizations for 3d on, when not debugging
+#		define OPT3D 1
+// COMPRT
+//#		define OPT3D 0
+#	endif // FSGR_STRICT_DEBUG==1
+#endif
+
+// enable/disable fine grid compression for finest level
+#if LBMDIM==3
+#define COMPRESSGRIDS 1
+#else 
+#define COMPRESSGRIDS 0
+#endif 
+
+//! threshold for level set fluid generation/isosurface
+#define LS_FLUIDTHRESHOLD 0.5
+
+//! invalid mass value for unused mass data
+#define MASS_INVALID -1000.0
+
+// empty/fill cells without fluid/empty NB's by inserting them into the full/empty lists?
+#define FSGR_LISTTRICK          1
+#define FSGR_LISTTTHRESHEMPTY   0.10
+#define FSGR_LISTTTHRESHFULL    0.90
+#define FSGR_MAGICNR            0.025
+//0.04
+
+//! maxmimum no. of grid levels
+#define FSGR_MAXNOOFLEVELS 5
+
+// helper for comparing floats with epsilon
+#define GFX_FLOATNEQ(x,y) ( ABS((x)-(y)) > (VECTOR_EPSILON) )
+#define LBM_FLOATNEQ(x,y) ( ABS((x)-(y)) > (10.0*LBM_EPSILON) )
+
+
+// macros for loops over all DFs
+#define FORDF0 for(int l= 0; l< LBM_DFNUM; ++l)
+#define FORDF1 for(int l= 1; l< LBM_DFNUM; ++l)
+// and with different loop var to prevent shadowing
+#define FORDF0M for(int m= 0; m< LBM_DFNUM; ++m)
+#define FORDF1M for(int m= 1; m< LBM_DFNUM; ++m)
+
+// aux. field indices (same for 2d)
+#define dFfrac 19
+#define dMass 20
+#define dFlux 21
+// max. no. of cell values for 3d
+#define dTotalNum 22
+
+// iso value defines
+// border for marching cubes
+#define ISOCORR 3
+
+
+// sirdude fix for solaris
+#if !defined(linux) && (defined (__sparc) || defined (__sparc__))
+#include <ieeefp.h>
+#ifndef expf
+#define expf exp
+#endif
+#endif
+
+
+/*****************************************************************************/
+/*! cell access classes */
+template<typename D>
+class UniformFsgrCellIdentifier : 
+	public CellIdentifierInterface 
+{
+	public:
+		//! which grid level?
+		int level;
+		//! location in grid
+		int x,y,z;
+
+		//! reset constructor
+		UniformFsgrCellIdentifier() :
+			x(0), y(0), z(0) { };
+
+		// implement CellIdentifierInterface
+		virtual string getAsString() {
+			std::ostringstream ret;
+			ret <<"{ i"<<x<<",j"<<y;
+			if(D::cDimension>2) ret<<",k"<<z;
+			ret <<" }";
+			return ret.str();
+		}
+
+		virtual bool equal(CellIdentifierInterface* other) {
+			//UniformFsgrCellIdentifier<D> *cid = dynamic_cast<UniformFsgrCellIdentifier<D> *>( other );
+			UniformFsgrCellIdentifier<D> *cid = (UniformFsgrCellIdentifier<D> *)( other );
+			if(!cid) return false;
+			if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true;
+			return false;
+		}
+};
+
+//! information needed for each level in the simulation
+class FsgrLevelData {
+public:
+	int id; // level number
+
+	//! node size on this level (geometric, in world coordinates, not simulation units!) 
+	LbmFloat nodeSize;
+	//! node size on this level in simulation units 
+	LbmFloat simCellSize;
+	//! quadtree node relaxation parameter 
+	LbmFloat omega;
+	//! size this level was advanced to 
+	LbmFloat time;
+	//! size of a single lbm step in time units on this level 
+	LbmFloat stepsize;
+	//! step count
+	int lsteps;
+	//! gravity force for this level
+	LbmVec gravity;
+	//! level array 
+	LbmFloat *mprsCells[2];
+	CellFlagType *mprsFlags[2];
+
+	//! smago params and precalculated values
+	LbmFloat lcsmago;
+	LbmFloat lcsmago_sqr;
+	LbmFloat lcnu;
+
+	// LES statistics per level
+	double avgOmega;
+	double avgOmegaCnt;
+
+	//! current set of dist funcs 
+	int setCurr;
+	//! target/other set of dist funcs 
+	int setOther;
+
+	//! mass&volume for this level
+	LbmFloat lmass;
+	LbmFloat lvolume;
+	LbmFloat lcellfactor;
+
+	//! local storage of mSizes
+	int lSizex, lSizey, lSizez;
+	int lOffsx, lOffsy, lOffsz;
+
+};
+
+
+
+/*****************************************************************************/
+/*! class for solving a LBM problem */
+template<class D>
+class LbmFsgrSolver : 
+	public D // this means, the solver is a lbmData object and implements the lbmInterface
+{
+
+	public:
+		//! Constructor 
+		LbmFsgrSolver();
+		//! Destructor 
+		virtual ~LbmFsgrSolver();
+		//! id string of solver
+		virtual string getIdString() { return string("FsgrSolver[") + D::getIdString(); }
+
+		//! initilize variables fom attribute list 
+		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 initialize( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ );
+
+#if LBM_USE_GUI==1
+		//! show simulation info (implement LbmSolverInterface pure virtual func)
+		virtual void debugDisplay(fluidDispSettings *set);
+#endif
+		
+		
+		// implement CellIterator<UniformFsgrCellIdentifier> interface
+		typedef UniformFsgrCellIdentifier<typename D::LbmCellContents> stdCellId;
+		virtual CellIdentifierInterface* getFirstCell( );
+		virtual void advanceCell( CellIdentifierInterface* );
+		virtual bool noEndCell( CellIdentifierInterface* );
+		virtual void deleteCellIterator( CellIdentifierInterface** );
+		virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos );
+		virtual int        getCellSet      ( CellIdentifierInterface* );
+		virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* );
+		virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* );
+		virtual int        getCellLevel    ( CellIdentifierInterface* );
+		virtual LbmFloat   getCellDensity  ( CellIdentifierInterface* ,int set);
+		virtual LbmVec     getCellVelocity ( CellIdentifierInterface* ,int set);
+		virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir);
+		virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set);
+		virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set);
+		virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set);
+		virtual LbmFloat   getEquilDf      ( int );
+		virtual int        getDfNum        ( );
+		// convert pointers
+		stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid);
+
+		//! perform geometry init (if switched on) 
+		bool initGeometryFlags();
+		//! init part for all freesurface testcases 
+		void initFreeSurfaces();
+		//! init density gradient if enabled
+		void initStandingFluidGradient();
+
+ 		/*! init a given cell with flag, density, mass and equilibrium dist. funcs */
+		inline void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
+		inline void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
+		inline void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
+
+		/*! perform a single LBM step */
+		virtual void step() { stepMain(); }
+		void stepMain();
+		void fineAdvance();
+		void coarseAdvance(int lev);
+		void coarseCalculateFluxareas(int lev);
+		// coarsen a given level (returns true if sth. was changed)
+		bool performRefinement(int lev);
+		bool performCoarsening(int lev);
+		//void oarseInterpolateToFineSpaceTime(int lev,LbmFloat t);
+		void interpolateFineFromCoarse(int lev,LbmFloat t);
+		void coarseRestrictFromFine(int lev);
+
+		/*! init particle positions */
+		virtual int initParticles(ParticleTracer *partt);
+		/*! move all particles */
+		virtual void advanceParticles(ParticleTracer *partt );
+
+
+		/*! debug object display (used e.g. for preview surface) */
+		virtual vector<ntlGeometryObject*> getDebugObjects();
+	
+		// gui/output debugging functions
+#if LBM_USE_GUI==1
+		virtual void debugDisplayNode(fluidDispSettings *dispset, CellIdentifierInterface* cell );
+		virtual void lbmDebugDisplay(fluidDispSettings *dispset);
+		virtual void lbmMarkedCellDisplay();
+#endif // LBM_USE_GUI==1
+		virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1);
+
+		//! for raytracing, preprocess
+		void prepareVisualization( void );
+
+		/*! type for cells */
+		typedef typename D::LbmCell LbmCell;
+		
+	protected:
+
+		//! internal quick print function (for debugging) 
+		void printLbmCell(int level, int i, int j, int k,int set);
+		// debugging use CellIterator interface to mark cell
+		void debugMarkCellCall(int level, int vi,int vj,int vk);
+
+		void mainLoop(int lev);
+		void adaptTimestep();
+		//! init mObjectSpeeds for current parametrization
+		void recalculateObjectSpeeds();
+		//! flag reinit step - always works on finest grid!
+		void reinitFlags( int workSet );
+		//! mass dist weights
+		LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l);
+		//! add point to mListNewInter list
+		inline void addToNewInterList( int ni, int nj, int nk );	
+		//! cell is interpolated from coarse level (inited into set, source sets are determined by t)
+		// inline, test?
+		void interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet,bool markNbs);
+
+		//! minimal and maximal z-coords (for 2D/3D loops)
+		int getForZMinBnd() { return 0; }
+		int getForZMin1()   { 
+			if(D::cDimension==2) return 0;
+			return 1; 
+		}
+
+		int getForZMaxBnd(int lev) { 
+			if(D::cDimension==2) return 1;
+			return mLevel[lev].lSizez -0;
+		}
+		int getForZMax1(int lev)   { 
+			if(D::cDimension==2) return 1;
+			return mLevel[lev].lSizez -1;
+		}
+
+
+		// member vars
+
+		//! mass calculated during streaming step
+		LbmFloat mCurrentMass;
+		LbmFloat mCurrentVolume;
+		LbmFloat mInitialMass;
+
+		//! count problematic cases, that occured so far...
+		int mNumProblems;
+
+		// average mlsups, count how many so far...
+		double mAvgMLSUPS;
+		double mAvgMLSUPSCnt;
+
+		//! Mcubes object for surface reconstruction 
+		IsoSurface *mpPreviewSurface;
+		float mSmoothSurface;
+		float mSmoothNormals;
+		
+		//! use time adaptivity? 
+		bool mTimeAdap;
+		//! domain boundary free/no slip type
+		string mDomainBound;
+		//! part slip value for domain
+		LbmFloat mDomainPartSlipValue;
+
+		//! output surface preview? if >0 yes, and use as reduzed size 
+		int mOutputSurfacePreview;
+		LbmFloat mPreviewFactor;
+		//! fluid vol height
+		LbmFloat mFVHeight;
+		LbmFloat mFVArea;
+		bool mUpdateFVHeight;
+
+		//! require some geo setup from the viz?
+		//int mGfxGeoSetup;
+		//! force quit for gfx
+		LbmFloat mGfxEndTime;
+		//! smoother surface initialization?
+		int mInitSurfaceSmoothing;
+
+		int mTimestepReduceLock;
+		int mTimeSwitchCounts;
+		//! total simulation time so far 
+		LbmFloat mSimulationTime;
+		//! smallest and largest step size so far 
+		LbmFloat mMinStepTime, mMaxStepTime;
+		//! track max. velocity
+		LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen;
+
+		//! list of the cells to empty at the end of the step 
+		vector<LbmPoint> mListEmpty;
+		//! list of the cells to make fluid at the end of the step 
+		vector<LbmPoint> mListFull;
+		//! list of new interface cells to init
+  	vector<LbmPoint> mListNewInter;
+		//! class for handling redist weights in reinit flag function
+		class lbmFloatSet {
+			public:
+				LbmFloat val[dTotalNum];
+				LbmFloat numNbs;
+		};
+		//! normalized vectors for all neighboring cell directions (for e.g. massdweight calc)
+		LbmVec mDvecNrm[27];
+		
+		
+		//! debugging
+		bool checkSymmetry(string idstring);
+		//! kepp track of max/min no. of filled cells
+		int mMaxNoCells, mMinNoCells;
+#ifndef USE_MSVC6FIXES
+		long long int mAvgNumUsedCells;
+#else
+		_int64 mAvgNumUsedCells;
+#endif
+
+		//! for interactive - how to drop drops?
+		int mDropMode;
+		LbmFloat mDropSize;
+		LbmVec mDropSpeed;
+		//! precalculated objects speeds for current parametrization
+		vector<LbmVec> mObjectSpeeds;
+		//! partslip bc. values for obstacle boundary conditions
+		vector<LbmFloat> mObjectPartslips;
+
+		//! get isofield weights
+		int mIsoWeightMethod;
+		float mIsoWeight[27];
+
+		// grid coarsening vars
+		
+		/*! vector for the data for each level */
+		FsgrLevelData mLevel[FSGR_MAXNOOFLEVELS];
+
+		/*! minimal and maximal refinement levels */
+		int mMaxRefine;
+
+		/*! df scale factors for level up/down */
+		LbmFloat mDfScaleUp, mDfScaleDown;
+
+		/*! precomputed cell area values */
+		LbmFloat mFsgrCellArea[27];
+
+		/*! LES C_smago paramter for finest grid */
+		float mInitialCsmago;
+		/*! LES stats for non OPT3D */
+		LbmFloat mDebugOmegaRet;
+
+		//! fluid stats
+		int mNumInterdCells;
+		int mNumInvIfCells;
+		int mNumInvIfTotal;
+		int mNumFsgrChanges;
+
+		//! debug function to disable standing f init
+		int mDisableStandingFluidInit;
+		//! debug function to force tadap syncing
+		int mForceTadapRefine;
+
+
+		// strict debug interface
+#		if FSGR_STRICT_DEBUG==1
+		int debLBMGI(int level, int ii,int ij,int ik, int is);
+		CellFlagType& debRFLAG(int level, int xx,int yy,int zz,int set);
+		CellFlagType& debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir);
+		CellFlagType& debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir);
+		int debLBMQI(int level, int ii,int ij,int ik, int is, int l);
+		LbmFloat& debQCELL(int level, int xx,int yy,int zz,int set,int l);
+		LbmFloat& debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l);
+		LbmFloat& debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l);
+		LbmFloat* debRACPNT(int level,  int ii,int ij,int ik, int is );
+		LbmFloat& debRAC(LbmFloat* s,int l);
+#		endif // FSGR_STRICT_DEBUG==1
+};
+
+
+
+/*****************************************************************************/
+// relaxation_macros
+
+
+
+// cell mark debugging
+#if FSGR_STRICT_DEBUG==10
+#define debugMarkCell(lev,x,y,z) \
+	errMsg("debugMarkCell",D::mName<<" step: "<<D::mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
+	debugMarkCellCall((lev),(x),(y),(z));
+#else // FSGR_STRICT_DEBUG==1
+#define debugMarkCell(lev,x,y,z) \
+	debugMarkCellCall((lev),(x),(y),(z));
+#endif // FSGR_STRICT_DEBUG==1
+
+
+// flag array defines -----------------------------------------------------------------------------------------------
+
+// lbm testsolver get index define
+#define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
+
+//! flag array acces macro
+#define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ]
+#define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set) ]
+#define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set) ]
+
+// array data layouts
+// standard array layout  -----------------------------------------------------------------------------------------------
+#define ALSTRING "Standard Array Layout"
+//#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
+#define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
+#define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
+#define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set, l)*dTotalNum +(l)])
+#define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set, l)*dTotalNum +(l)])
+
+#define QCELLSTEP dTotalNum
+#define _RACPNT(level, ii,ij,ik, is )  &QCELL(level,ii,ij,ik,is,0)
+#define _RAC(s,l) (s)[(l)]
+
+
+#if FSGR_STRICT_DEBUG==1
+
+#define LBMGI(level,ii,ij,ik, is)                 debLBMGI(level,ii,ij,ik, is)         
+#define RFLAG(level,xx,yy,zz,set)                 debRFLAG(level,xx,yy,zz,set)            
+#define RFLAG_NB(level,xx,yy,zz,set, dir)         debRFLAG_NB(level,xx,yy,zz,set, dir)    
+#define RFLAG_NBINV(level,xx,yy,zz,set, dir)      debRFLAG_NBINV(level,xx,yy,zz,set, dir) 
+
+#define LBMQI(level,ii,ij,ik, is, l)              debLBMQI(level,ii,ij,ik, is, l)         
+#define QCELL(level,xx,yy,zz,set,l)               debQCELL(level,xx,yy,zz,set,l)         
+#define QCELL_NB(level,xx,yy,zz,set, dir,l)       debQCELL_NB(level,xx,yy,zz,set, dir,l)
+#define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    debQCELL_NBINV(level,xx,yy,zz,set, dir,l)
+#define RACPNT(level, ii,ij,ik, is )              debRACPNT(level, ii,ij,ik, is )          
+#define RAC(s,l)                            			debRAC(s,l)                  
+
+#else // FSGR_STRICT_DEBUG==1
+
+#define LBMGI(level,ii,ij,ik, is)                 _LBMGI(level,ii,ij,ik, is)         
+#define RFLAG(level,xx,yy,zz,set)                 _RFLAG(level,xx,yy,zz,set)            
+#define RFLAG_NB(level,xx,yy,zz,set, dir)         _RFLAG_NB(level,xx,yy,zz,set, dir)    
+#define RFLAG_NBINV(level,xx,yy,zz,set, dir)      _RFLAG_NBINV(level,xx,yy,zz,set, dir) 
+
+#define LBMQI(level,ii,ij,ik, is, l)              _LBMQI(level,ii,ij,ik, is, l)         
+#define QCELL(level,xx,yy,zz,set,l)               _QCELL(level,xx,yy,zz,set,l)         
+#define QCELL_NB(level,xx,yy,zz,set, dir,l)       _QCELL_NB(level,xx,yy,zz,set, dir, l)
+#define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    _QCELL_NBINV(level,xx,yy,zz,set, dir,l)
+#define RACPNT(level, ii,ij,ik, is )              _RACPNT(level, ii,ij,ik, is )          
+#define RAC(s,l)                                  _RAC(s,l)                  
+
+#endif // FSGR_STRICT_DEBUG==1
+
+// general defines -----------------------------------------------------------------------------------------------
+
+#define TESTFLAG(flag, compflag) ((flag & compflag)==compflag)
+
+#if LBMDIM==2
+#define dC 0
+#define dN 1
+#define dS 2
+#define dE 3
+#define dW 4
+#define dNE 5
+#define dNW 6
+#define dSE 7
+#define dSW 8
+#define LBM_DFNUM 9
+#else
+// direction indices
+#define dC 0
+#define dN 1
+#define dS 2
+#define dE 3
+#define dW 4
+#define dT 5
+#define dB 6
+#define dNE 7
+#define dNW 8
+#define dSE 9
+#define dSW 10
+#define dNT 11
+#define dNB 12
+#define dST 13
+#define dSB 14
+#define dET 15
+#define dEB 16
+#define dWT 17
+#define dWB 18
+#define LBM_DFNUM 19
+#endif
+//? #define dWB 18
+
+// default init for dFlux values
+#define FLUX_INIT 0.5f * (float)(D::cDfNum)
+
+// only for non DF dir handling!
+#define dNET 19
+#define dNWT 20
+#define dSET 21
+#define dSWT 22
+#define dNEB 23
+#define dNWB 24
+#define dSEB 25
+#define dSWB 26
+
+//! fill value for boundary cells
+#define BND_FILL 0.0
+
+#define DFL1 (1.0/ 3.0)
+#define DFL2 (1.0/18.0)
+#define DFL3 (1.0/36.0)
+
+// loops over _all_ cells (including boundary layer)
+#define FSGR_FORIJK_BOUNDS(leveli) \
+	for(int k= getForZMinBnd(); k< getForZMaxBnd(leveli); ++k) \
+   for(int j=0;j<mLevel[leveli].lSizey-0;++j) \
+    for(int i=0;i<mLevel[leveli].lSizex-0;++i) \
+	
+// loops over _only inner_ cells 
+#define FSGR_FORIJK1(leveli) \
+	for(int k= getForZMin1(); k< getForZMax1(leveli); ++k) \
+   for(int j=1;j<mLevel[leveli].lSizey-1;++j) \
+    for(int i=1;i<mLevel[leveli].lSizex-1;++i) \
+
+// relaxation_macros end
+
+
+/*****************************************************************************/
+/* init a given cell with flag, density, mass and equilibrium dist. funcs */
+
+template<class D>
+void LbmFsgrSolver<D>::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
+	CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask;
+	RFLAG(level,xx,yy,zz,set) = newflag | pers;
+}
+
+template<class D>
+void 
+LbmFsgrSolver<D>::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
+  /* init eq. dist funcs */
+	LbmFloat *ecel;
+	int workSet = mLevel[level].setCurr;
+
+	ecel = RACPNT(level, i,j,k, workSet);
+	FORDF0 { RAC(ecel, l) = D::dfEquil[l] * rho; }
+	RAC(ecel, dMass) = mass;
+	RAC(ecel, dFfrac) = mass/rho;
+	RAC(ecel, dFlux) = FLUX_INIT;
+	//RFLAG(level, i,j,k, workSet)= flag;
+	changeFlag(level, i,j,k, workSet, flag);
+
+  workSet ^= 1;
+	changeFlag(level, i,j,k, workSet, flag);
+	return;
+}
+
+template<class D>
+void 
+LbmFsgrSolver<D>::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
+	LbmFloat *ecel;
+	int workSet = mLevel[level].setCurr;
+
+	ecel = RACPNT(level, i,j,k, workSet);
+	FORDF0 { RAC(ecel, l) = D::getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
+	RAC(ecel, dMass) = mass;
+	RAC(ecel, dFfrac) = mass/rho;
+	RAC(ecel, dFlux) = FLUX_INIT;
+	//RFLAG(level, i,j,k, workSet) = flag;
+	changeFlag(level, i,j,k, workSet, flag);
+
+  workSet ^= 1;
+	changeFlag(level, i,j,k, workSet, flag);
+	return;
+}
+
+#define LBMFSGRSOLVER_H
+#endif
+
+

intern/elbeem/intern/solver_dimenions.h

@@ -0,0 +1,21 @@
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * All code distributed as part of El'Beem is covered by the version 2 of the 
+ * GNU General Public License. See the file COPYING for details.
+ * Copyright 2003-2005 Nils Thuerey
+ *
+ * Combined 2D/3D Lattice Boltzmann Solver auxiliary classes
+ * 
+ *****************************************************************************/
+#ifndef LBMHEADER_H
+
+/* LBM Files */
+#include "solver_interface.h"
+
+
+#define LBMHEADER_H
+#endif
+
+
+

intern/elbeem/intern/solver_interface.cpp

@@ -11,12 +11,9 @@
  *****************************************************************************/
 
 /* LBM Files */ 
-#include "lbmdimensions.h" 
+#include "solver_interface.h" 
-#include "lbminterface.h" 
-#include "lbmfunctions.h" 
 #include "ntl_scene.h"
 #include "ntl_ray.h"
-#include "typeslbm.h"
 
 
 /*****************************************************************************/
@@ -401,6 +398,7 @@ void LbmSolverInterface::initGeoTree(int id) {
 	mpGiObjects = scene->getObjects();
 	mGiObjInside.resize( mpGiObjects->size() );
 	mGiObjDistance.resize( mpGiObjects->size() );
+	mGiObjSecondDist.resize( mpGiObjects->size() );
 	for(size_t i=0; i<mpGiObjects->size(); i++) { 
 		if((*mpGiObjects)[i]->getGeoInitIntersect()) mAccurateGeoinit=true;
 	}
@@ -436,6 +434,7 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int
 	for(size_t i=0; i<mGiObjInside.size(); i++) { 
 		mGiObjInside[i] = 0; 
 		mGiObjDistance[i] = -1.0; 
+		mGiObjSecondDist[i] = -1.0; 
 	}
 	// if not inside, return distance to first hit
 	gfxReal firstHit=-1.0;
@@ -517,6 +516,142 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int
 		return false;
 	}
 }
+bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, const gfxReal halfCellsize, bool &thinHit, bool recurse) {
+	// shift ve ctors to avoid rounding errors
+	org += ntlVec3Gfx(0.0001); //?
+	OId = -1;
+	ntlRay ray(org, dir, 0, 1.0, mpGlob);
+	//int insCnt = 0;
+	bool done = false;
+	bool inside = false;
+	for(size_t i=0; i<mGiObjInside.size(); i++) { 
+		mGiObjInside[i] = 0; 
+		mGiObjDistance[i] = -1.0; 
+		mGiObjSecondDist[i] = -1.0; 
+	}
+	// if not inside, return distance to first hit
+	gfxReal firstHit=-1.0;
+	int     firstOId = -1;
+	thinHit = false;
+	
+	if(mAccurateGeoinit) {
+		while(!done) {
+			// find first inside intersection
+			ntlTriangle *triIns = NULL;
+			distance = -1.0;
+			ntlVec3Gfx normal(0.0);
+			mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
+			if(triIns) {
+				ntlVec3Gfx norg = ray.getOrigin() + ray.getDirection()*distance;
+				LbmFloat orientation = dot(normal, dir);
+				OId = triIns->getObjectId();
+				if(orientation<=0.0) {
+					// outside hit
+					normal *= -1.0;
+					//mGiObjDistance[OId] = -1.0;
+					//errMsg("IIO"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
+				} else {
+					// inside hit
+					//if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
+					//errMsg("III"," oid:"<<OId<<" org"<<org<<" norg"<<norg);
+					if(mGiObjInside[OId]==1) {
+						// second inside hit
+						if(mGiObjSecondDist[OId]<0.0) mGiObjSecondDist[OId] = distance;
+					}
+				}
+				mGiObjInside[OId]++;
+				// always store first hit for thin obj init
+				if(mGiObjDistance[OId]<0.0) mGiObjDistance[OId] = distance;
+
+				norg += normal * getVecEpsilon();
+				ray = ntlRay(norg, dir, 0, 1.0, mpGlob);
+				// remember first hit distance, in case we're not 
+				// inside anything
+				if(firstHit<0.0) {
+					firstHit = distance;
+					firstOId = OId;
+				}
+			} else {
+				// no more intersections... return false
+				done = true;
+				//if(insCnt%2) inside=true;
+			}
+		}
+
+		distance = -1.0;
+		// standard inside check
+		for(size_t i=0; i<mGiObjInside.size(); i++) {
+			if(((mGiObjInside[i]%2)==1)&&(mGiObjDistance[i]>0.0)) {
+				if(  (distance<0.0)                             || // first intersection -> good
+					  ((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
+						) {						
+					distance = mGiObjDistance[i];
+					OId = i;
+					inside = true;
+				} 
+			}
+		}
+		// now check for thin hits
+		if(!inside) {
+			distance = -1.0;
+			for(size_t i=0; i<mGiObjInside.size(); i++) {
+				if((mGiObjInside[i]>=2)&&(mGiObjDistance[i]>0.0)&&(mGiObjSecondDist[i]>0.0)&&
+					 (mGiObjDistance[i]<1.0*halfCellsize)&&(mGiObjSecondDist[i]<2.0*halfCellsize) ) {
+					if(  (distance<0.0)                             || // first intersection -> good
+							((distance>0.0)&&(distance>mGiObjDistance[i])) // more than one intersection -> use closest one
+							) {						
+						distance = mGiObjDistance[i];
+						OId = i;
+						inside = true;
+						thinHit = true;
+					} 
+				}
+			}
+		}
+		if(!inside) {
+			// check for hit in this cell, opposite to current dir (only recurse once)
+			if(recurse) {
+				gfxReal r_distance;
+				int r_OId;
+				bool ret = geoInitCheckPointInside(org, dir*-1.0, flags, r_OId, r_distance, halfCellsize, thinHit, false);
+				if((ret)&&(thinHit)) {
+					OId = r_OId;
+					distance = 0.0; 
+					return true;
+				}
+			}
+		}
+		// really no hit...
+		if(!inside) {
+			distance = firstHit;
+			OId = firstOId;
+			/*if((mGiObjDistance[OId]>0.0)&&(mGiObjSecondDist[OId]>0.0)) {
+				const gfxReal thisdist = mGiObjSecondDist[OId]-mGiObjDistance[OId];
+				// dont walk over this cell...
+				if(thisdist<halfCellsize) distance-=2.0*halfCellsize;
+			} // ? */
+		}
+		//errMsg("CHIII","i"<<inside<<"  fh"<<firstHit<<" fo"<<firstOId<<" - h"<<distance<<" o"<<OId);
+
+		return inside;
+	} else {
+
+		// find first inside intersection
+		ntlTriangle *triIns = NULL;
+		distance = -1.0;
+		ntlVec3Gfx normal(0.0);
+		mpGiTree->intersect(ray,distance,normal, triIns, flags, true);
+		if(triIns) {
+			// check outside intersect
+			LbmFloat orientation = dot(normal, dir);
+			if(orientation<=0.0) return false;
+
+			OId = triIns->getObjectId();
+			return true;
+		}
+		return false;
+	}
+}
 
 /*****************************************************************************/
 /*! get max. velocity of all objects to initialize as fluid regions */

intern/elbeem/intern/solver_interface.h

@@ -288,6 +288,8 @@ class LbmSolverInterface
 		void freeGeoTree();
 		/*! check for a certain flag type at position org (needed for e.g. quadtree refinement) */
 		bool geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance);
+		bool geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, 
+				const gfxReal halfCellsize, bool &thinHit, bool recurse);
 		/*! set render globals, for scene/tree access */
 		void setRenderGlobals(ntlRenderGlobals *glob) { mpGlob = glob; };
 		/*! get max. velocity of all objects to initialize as fluid regions */
@@ -378,6 +380,13 @@ class LbmSolverInterface
 		virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set) = 0;
 		virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set) = 0;
 
+		// gui/output debugging functions
+#if LBM_USE_GUI==1
+		virtual void debugDisplayNode(fluidDispSettings *dispset, CellIdentifier cell ) = 0;
+		virtual void lbmDebugDisplay(fluidDispSettings *dispset) = 0;
+		virtual void lbmMarkedCellDisplay() = 0;
+#endif // LBM_USE_GUI==1
+		virtual void debugPrintNodeInfo(CellIdentifier cell, int forceSet=-1) = 0;
 
 		// debugging cell marker functions
 
@@ -500,6 +509,7 @@ class LbmSolverInterface
 		vector<int> mGiObjInside;
 		/*! inside which objects? */
 		vector<gfxReal> mGiObjDistance;
+		vector<gfxReal> mGiObjSecondDist;
 		/*! remember globals */
 		ntlRenderGlobals *mpGlob;
 
@@ -509,6 +519,375 @@ class LbmSolverInterface
 };
 
 
+//! shorten static const definitions
+#define STCON static const
+
+
+/*****************************************************************************/
+/*! class for solver templating - 3D implementation */
+class LbmD3Q19 {
+
+	public:
+
+		// constructor, init interface
+		LbmD3Q19() {};
+		// virtual destructor 
+		virtual ~LbmD3Q19() {};
+		//! id string of solver
+		string getIdString() { return string("3D"); }
+
+		//! how many dimensions?
+		STCON int cDimension;
+
+		// Wi factors for collide step 
+		STCON LbmFloat cCollenZero;
+		STCON LbmFloat cCollenOne;
+		STCON LbmFloat cCollenSqrtTwo;
+
+		//! threshold value for filled/emptied cells 
+		STCON LbmFloat cMagicNr2;
+		STCON LbmFloat cMagicNr2Neg;
+		STCON LbmFloat cMagicNr;
+		STCON LbmFloat cMagicNrNeg;
+
+		//! size of a single set of distribution functions 
+		STCON int    cDfNum;
+		//! direction vector contain vecs for all spatial dirs, even if not used for LBM model
+		STCON int    cDirNum;
+
+		//! distribution functions directions 
+		typedef enum {
+			 cDirInv=  -1,
+			 cDirC  =  0,
+			 cDirN  =  1,
+			 cDirS  =  2,
+			 cDirE  =  3,
+			 cDirW  =  4,
+			 cDirT  =  5,
+			 cDirB  =  6,
+			 cDirNE =  7,
+			 cDirNW =  8,
+			 cDirSE =  9,
+			 cDirSW = 10,
+			 cDirNT = 11,
+			 cDirNB = 12,
+			 cDirST = 13,
+			 cDirSB = 14,
+			 cDirET = 15,
+			 cDirEB = 16,
+			 cDirWT = 17,
+			 cDirWB = 18
+		} dfDir;
+
+		/* Vector Order 3D:
+		 *  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
+		 *  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
+		 *  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
+		 *  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1
+		 */
+
+		/*! name of the dist. function 
+			 only for nicer output */
+		STCON char* dfString[ 19 ];
+
+		/*! index of normal dist func, not used so far?... */
+		STCON int dfNorm[ 19 ];
+
+		/*! index of inverse dist func, not fast, but useful... */
+		STCON int dfInv[ 19 ];
+
+		/*! index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefX[ 19 ];
+		/*! index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefY[ 19 ];
+		/*! index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefZ[ 19 ];
+
+		/*! dist func vectors */
+		STCON int dfVecX[ 27 ];
+		STCON int dfVecY[ 27 ];
+		STCON int dfVecZ[ 27 ];
+
+		/*! arrays as before with doubles */
+		STCON LbmFloat dfDvecX[ 27 ];
+		STCON LbmFloat dfDvecY[ 27 ];
+		STCON LbmFloat dfDvecZ[ 27 ];
+
+		/*! principal directions */
+		STCON int princDirX[ 2*3 ];
+		STCON int princDirY[ 2*3 ];
+		STCON int princDirZ[ 2*3 ];
+
+		/*! vector lengths */
+		STCON LbmFloat dfLength[ 19 ];
+
+		/*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
+		static LbmFloat dfEquil[ 19 ];
+
+		/*! arrays for les model coefficients */
+		static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
+		static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
+
+}; // LbmData3D
+
+
+
+/*****************************************************************************/
+//! class for solver templating - 2D implementation 
+class LbmD2Q9 {
+	
+	public:
+
+		// constructor, init interface
+		LbmD2Q9() {};
+		// virtual destructor 
+		virtual ~LbmD2Q9() {};
+		//! id string of solver
+		string getIdString() { return string("2D"); }
+
+		//! how many dimensions?
+		STCON int cDimension;
+
+		//! Wi factors for collide step 
+		STCON LbmFloat cCollenZero;
+		STCON LbmFloat cCollenOne;
+		STCON LbmFloat cCollenSqrtTwo;
+
+		//! threshold value for filled/emptied cells 
+		STCON LbmFloat cMagicNr2;
+		STCON LbmFloat cMagicNr2Neg;
+		STCON LbmFloat cMagicNr;
+		STCON LbmFloat cMagicNrNeg;
+
+		//! size of a single set of distribution functions 
+		STCON int    cDfNum;
+		STCON int    cDirNum;
+
+		//! distribution functions directions 
+		typedef enum {
+			 cDirInv=  -1,
+			 cDirC  =  0,
+			 cDirN  =  1,
+			 cDirS  =  2,
+			 cDirE  =  3,
+			 cDirW  =  4,
+			 cDirNE =  5,
+			 cDirNW =  6,
+			 cDirSE =  7,
+			 cDirSW =  8
+		} dfDir;
+
+		/* Vector Order 2D:
+		 * 0  1 2  3  4  5  6 7  8
+		 * 0, 0,0, 1,-1, 1,-1,1,-1 
+		 * 0, 1,-1, 0,0, 1,1,-1,-1  */
+
+		/* name of the dist. function 
+			 only for nicer output */
+		STCON char* dfString[ 9 ];
+
+		/* index of normal dist func, not used so far?... */
+		STCON int dfNorm[ 9 ];
+
+		/* index of inverse dist func, not fast, but useful... */
+		STCON int dfInv[ 9 ];
+
+		/* index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefX[ 9 ];
+		/* index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefY[ 9 ];
+		/* index of x reflected dist func for free slip, not valid for all DFs... */
+		STCON int dfRefZ[ 9 ];
+
+		/* dist func vectors */
+		STCON int dfVecX[ 9 ];
+		STCON int dfVecY[ 9 ];
+		/* Z, 2D values are all 0! */
+		STCON int dfVecZ[ 9 ];
+
+		/* arrays as before with doubles */
+		STCON LbmFloat dfDvecX[ 9 ];
+		STCON LbmFloat dfDvecY[ 9 ];
+		/* Z, 2D values are all 0! */
+		STCON LbmFloat dfDvecZ[ 9 ];
+
+		/*! principal directions */
+		STCON int princDirX[ 2*2 ];
+		STCON int princDirY[ 2*2 ];
+		STCON int princDirZ[ 2*2 ];
+
+		/* vector lengths */
+		STCON LbmFloat dfLength[ 9 ];
+
+		/* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
+		static LbmFloat dfEquil[ 9 ];
+
+		/*! arrays for les model coefficients */
+		static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
+		static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
+
+}; // LbmData3D
+
+
+
+// lbmdimensions
+
+// not needed hereafter
+#undef STCON
+
+
+
+/*****************************************************************************/
+//! class for solver templating - lbgk (srt) model implementation 
+template<class DQ>
+class LbmModelLBGK : public DQ , public LbmSolverInterface {
+	public:
+
+		/*! type for cells contents, needed for cell id interface */
+		typedef DQ LbmCellContents;
+		/*! type for cells */
+		typedef LbmCellTemplate< LbmCellContents > LbmCell;
+
+		// constructor
+		LbmModelLBGK() : DQ(), LbmSolverInterface() {};
+		// virtual destructor 
+		virtual ~LbmModelLBGK() {};
+		//! id string of solver
+		string getIdString() { return DQ::getIdString() + string("lbgk]"); }
+
+		/*! calculate length of velocity vector */
+		static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) {
+			return ((ux)*DQ::dfDvecX[l]+(uy)*DQ::dfDvecY[l]+(uz)*DQ::dfDvecZ[l]);
+		};
+
+		/*! calculate equilibrium DF for given values */
+		static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz) {
+			LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 
+			return( DQ::dfLength[l] *( 
+						+ rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
+						+ 3.0 *tmp 
+						+ 9.0/2.0 *(tmp*tmp) ) 
+					);
+		};
+
+
+		// input mux etc. as acceleration
+		// outputs rho,ux,uy,uz
+		/*inline void collideArrays_org(LbmFloat df[19], 				
+				LbmFloat &outrho, // out only!
+				// velocity modifiers (returns actual velocity!)
+				LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
+				LbmFloat omega
+			) {
+			LbmFloat rho=df[0]; 
+			LbmFloat ux = mux;
+			LbmFloat uy = muy;
+			LbmFloat uz = muz;
+			for(int l=1; l<DQ::cDfNum; l++) { 
+				rho += df[l]; 
+				ux  += (DQ::dfDvecX[l]*df[l]); 
+				uy  += (DQ::dfDvecY[l]*df[l]);  
+				uz  += (DQ::dfDvecZ[l]*df[l]);  
+			}  
+			for(int l=0; l<DQ::cDfNum; l++) { 
+				//LbmFloat tmp = (ux*DQ::dfDvecX[l]+uy*DQ::dfDvecY[l]+uz*DQ::dfDvecZ[l]); 
+				df[l] = (1.0-omega ) * df[l] + omega * ( getCollideEq(l,rho,ux,uy,uz) ); 
+			}  
+
+			mux = ux;
+			muy = uy;
+			muz = uz;
+			outrho = rho;
+		};*/
+		
+		// LES functions
+		inline LbmFloat getLesNoneqTensorCoeff(
+				LbmFloat df[], 				
+				LbmFloat feq[] ) {
+			LbmFloat Qo = 0.0;
+			for(int m=0; m< ((DQ::cDimension*DQ::cDimension)-DQ::cDimension)/2 ; m++) { 
+				LbmFloat qadd = 0.0;
+				for(int l=1; l<DQ::cDfNum; l++) { 
+					if(DQ::lesCoeffOffdiag[m][l]==0.0) continue;
+					qadd += DQ::lesCoeffOffdiag[m][l]*(df[l]-feq[l]);
+				}
+				Qo += (qadd*qadd);
+			}
+			Qo *= 2.0; // off diag twice
+			for(int m=0; m<DQ::cDimension; m++) { 
+				LbmFloat qadd = 0.0;
+				for(int l=1; l<DQ::cDfNum; l++) { 
+					if(DQ::lesCoeffDiag[m][l]==0.0) continue;
+					qadd += DQ::lesCoeffDiag[m][l]*(df[l]-feq[l]);
+				}
+				Qo += (qadd*qadd);
+			}
+			Qo = sqrt(Qo);
+			return Qo;
+		}
+		inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo) {
+			const LbmFloat tau = 1.0/omega;
+			const LbmFloat nu = (2.0*tau-1.0) * (1.0/6.0);
+			const LbmFloat C = csmago;
+			const LbmFloat Csqr = C*C;
+			LbmFloat S = -nu + sqrt( nu*nu + 18.0*Csqr*Qo ) / (6.0*Csqr);
+			return( 1.0/( 3.0*( nu+Csqr*S ) +0.5 ) );
+		}
+
+		// "normal" collision
+		inline void collideArrays(LbmFloat df[], 				
+				LbmFloat &outrho, // out only!
+				// velocity modifiers (returns actual velocity!)
+				LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
+				LbmFloat omega, LbmFloat csmago, LbmFloat *newOmegaRet = NULL
+			) {
+			LbmFloat rho=df[0]; 
+			LbmFloat ux = mux;
+			LbmFloat uy = muy;
+			LbmFloat uz = muz; 
+			for(int l=1; l<DQ::cDfNum; l++) { 
+				rho += df[l]; 
+				ux  += (DQ::dfDvecX[l]*df[l]); 
+				uy  += (DQ::dfDvecY[l]*df[l]);  
+				uz  += (DQ::dfDvecZ[l]*df[l]);  
+			}  
+			LbmFloat feq[19];
+			for(int l=0; l<DQ::cDfNum; l++) { 
+				feq[l] = getCollideEq(l,rho,ux,uy,uz); 
+			}
+
+			LbmFloat omegaNew;
+			if(csmago>0.0) {
+				LbmFloat Qo = getLesNoneqTensorCoeff(df,feq);
+				omegaNew = getLesOmega(omega,csmago,Qo);
+			} else {
+				omegaNew = omega; // smago off...
+			}
+			if(newOmegaRet) *newOmegaRet=omegaNew; // return value for stats
+
+			for(int l=0; l<DQ::cDfNum; l++) { 
+				df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
+			}  
+
+			mux = ux;
+			muy = uy;
+			muz = uz;
+			outrho = rho;
+		};
+
+}; // LBGK
+
+#ifdef LBMMODEL_DEFINED
+// force compiler error!
+ERROR - Dont include several LBM models at once...
+#endif
+#define LBMMODEL_DEFINED 1
+
+
+typedef LbmModelLBGK<  LbmD2Q9 > LbmBGK2D;
+typedef LbmModelLBGK< LbmD3Q19 > LbmBGK3D;
+
+
 //! helper function to convert flag to string (for debuggin)
 string convertCellFlagType2String( CellFlagType flag );
 string convertSingleFlag2String(CellFlagType cflag);

intern/elbeem/intern/solver_util.cpp

@@ -0,0 +1,861 @@
+/******************************************************************************
+ *
+ * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
+ * Copyright 2003,2004,2005 Nils Thuerey
+ *
+ * Standard LBM Factory implementation
+ *
+ *****************************************************************************/
+
+#include "solver_class.h"
+
+//! for raytracing
+template<class D>
+void LbmFsgrSolver<D>::prepareVisualization( void ) {
+	int lev = mMaxRefine;
+	int workSet = mLevel[lev].setCurr;
+
+	//make same prepareVisualization and getIsoSurface...
+#if LBMDIM==2
+	// 2d, place in the middle of isofield slice (k=2)
+#  define ZKD1 0
+	// 2d z offset = 2, lbmGetData adds 1, so use one here
+#  define ZKOFF 1
+	// reset all values...
+	for(int k= 0; k< 5; ++k) 
+   for(int j=0;j<mLevel[lev].lSizey-0;j++) 
+    for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+		*D::mpIso->lbmGetData(i,j,ZKOFF)=0.0;
+	}
+#else // LBMDIM==2
+	// 3d, use normal bounds
+#  define ZKD1 1
+#  define ZKOFF k
+	// reset all values...
+	for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
+   for(int j=0;j<mLevel[lev].lSizey-0;j++) 
+    for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+		*D::mpIso->lbmGetData(i,j,ZKOFF)=0.0;
+	}
+#endif // LBMDIM==2
+
+
+	
+	// add up...
+	float val = 0.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++) {
+
+		//continue; // OFF DEBUG
+		if(RFLAG(lev, i,j,k,workSet)&(CFBnd|CFEmpty)) {
+			continue;
+		} else
+		if( (RFLAG(lev, i,j,k,workSet)&CFInter) && (!(RFLAG(lev, i,j,k,workSet)&CFNoNbEmpty)) ){
+			// no empty nb interface cells are treated as full
+			val =  (QCELL(lev, i,j,k,workSet, dFfrac)); 
+			//if( (!(RFLAG(lev, i,j,k,workSet)&CFNoBndFluid)) &&(RFLAG(lev, i,j,k,workSet)&CFNoNbFluid)){ val += D::mIsoValue; }
+		} else {
+			// fluid?
+			val = 1.0; ///27.0;
+		} // */
+
+		*D::mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
+		*D::mpIso->lbmGetData( i   , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 
+		*D::mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); 
+										
+		*D::mpIso->lbmGetData( i-1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); 
+		*D::mpIso->lbmGetData( i   , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); 
+		*D::mpIso->lbmGetData( i+1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); 
+										
+		*D::mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); 
+		*D::mpIso->lbmGetData( i   , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); 
+		*D::mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); 
+										
+										
+		*D::mpIso->lbmGetData( i-1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[9] ); 
+		*D::mpIso->lbmGetData( i   , j-1  ,ZKOFF  ) += ( val * mIsoWeight[10] ); 
+		*D::mpIso->lbmGetData( i+1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[11] ); 
+																	
+		*D::mpIso->lbmGetData( i-1 , j    ,ZKOFF  ) += ( val * mIsoWeight[12] ); 
+		*D::mpIso->lbmGetData( i   , j    ,ZKOFF  ) += ( val * mIsoWeight[13] ); 
+		*D::mpIso->lbmGetData( i+1 , j    ,ZKOFF  ) += ( val * mIsoWeight[14] ); 
+																	
+		*D::mpIso->lbmGetData( i-1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[15] ); 
+		*D::mpIso->lbmGetData( i   , j+1  ,ZKOFF  ) += ( val * mIsoWeight[16] ); 
+		*D::mpIso->lbmGetData( i+1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[17] ); 
+										
+										
+		*D::mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); 
+		*D::mpIso->lbmGetData( i   , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); 
+		*D::mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); 
+																 
+		*D::mpIso->lbmGetData( i-1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); 
+		*D::mpIso->lbmGetData( i   , j   ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); 
+		*D::mpIso->lbmGetData( i+1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); 
+																 
+		*D::mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); 
+		*D::mpIso->lbmGetData( i   , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); 
+		*D::mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
+	}
+
+	// update preview, remove 2d?
+	if(mOutputSurfacePreview) {
+		//int previewSize = mOutputSurfacePreview;
+		int pvsx = (int)(mPreviewFactor*D::mSizex);
+		int pvsy = (int)(mPreviewFactor*D::mSizey);
+		int pvsz = (int)(mPreviewFactor*D::mSizez);
+		//float scale = (float)D::mSizex / previewSize;
+		LbmFloat scalex = (LbmFloat)D::mSizex/(LbmFloat)pvsx;
+		LbmFloat scaley = (LbmFloat)D::mSizey/(LbmFloat)pvsy;
+		LbmFloat scalez = (LbmFloat)D::mSizez/(LbmFloat)pvsz;
+		for(int k= 0; k< ((D::cDimension==3) ? (pvsz-1):1) ; ++k) 
+   		for(int j=0;j< pvsy;j++) 
+    		for(int i=0;i< pvsx;i++) {
+					*mpPreviewSurface->lbmGetData(i,j,k) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) );
+				}
+		// set borders again...
+		for(int k= 0; k< ((D::cDimension == 3) ? (pvsz-1):1) ; ++k) {
+			for(int j=0;j< pvsy;j++) {
+				*mpPreviewSurface->lbmGetData(0,j,k) = *D::mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) );
+				*mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *D::mpIso->lbmGetData( D::mSizex-1, (int)(j*scaley), (int)(k*scalez) );
+			}
+			for(int i=0;i< pvsx;i++) {
+				*mpPreviewSurface->lbmGetData(i,0,k) = *D::mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) );
+				*mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *D::mpIso->lbmGetData( (int)(i*scalex), D::mSizey-1, (int)(k*scalez) );
+			}
+		}
+		if(D::cDimension == 3) {
+			// only for 3D
+			for(int j=0;j<pvsy;j++)
+				for(int i=0;i<pvsx;i++) {      
+					*mpPreviewSurface->lbmGetData(i,j,0) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0);
+					*mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , D::mSizez-1);
+				}
+		} // borders done...
+	}
+
+	// correction
+	return;
+}
+
+
+
+
+/*****************************************************************************
+ * move the particles
+ * uses updated velocities from mSetOther
+ *****************************************************************************/
+template<class D>
+void LbmFsgrSolver<D>::advanceParticles(ParticleTracer *partt ) { 
+	partt = NULL; // remove warning
+}
+
+
+/******************************************************************************
+ * reset particle positions to default
+ *****************************************************************************/
+/*! init particle positions */
+template<class D>
+int LbmFsgrSolver<D>::initParticles(ParticleTracer *partt) { 
+	partt = NULL; // remove warning
+	return 0;
+}
+
+
+/*! init particle positions */
+template<class D>
+void LbmFsgrSolver<D>::recalculateObjectSpeeds() {
+	int numobjs = (int)(D::mpGiObjects->size());
+	// note - (numobjs + 1) is entry for domain settings
+	if(numobjs>255-1) {
+		errFatal("LbmFsgrSolver::recalculateObjectSpeeds","More than 256 objects currently not supported...",SIMWORLD_INITERROR);
+		return;
+	}
+	mObjectSpeeds.resize(numobjs+1);
+	for(int i=0; i<(int)(numobjs+0); i++) {
+		mObjectSpeeds[i] = vec2L(D::mpParam->calculateLattVelocityFromRw( vec2P( (*D::mpGiObjects)[i]->getInitialVelocity() )));
+		//errMsg("recalculateObjectSpeeds","id"<<i<<" set to "<< mObjectSpeeds[i]<<", unscaled:"<< (*D::mpGiObjects)[i]->getInitialVelocity() );
+	}
+
+	// also reinit part slip values here
+	mObjectPartslips.resize(numobjs+1);
+	for(int i=0; i<(int)(numobjs+0); i++) {
+		mObjectPartslips[i] = (LbmFloat)(*D::mpGiObjects)[i]->getGeoPartSlipValue();
+	}
+	//errMsg("GEOIN"," dm set "<<mDomainPartSlipValue);
+	mObjectPartslips[numobjs] = mDomainPartSlipValue;
+}
+
+/*****************************************************************************/
+/*! internal quick print function (for debugging) */
+/*****************************************************************************/
+template<class D>
+void 
+LbmFsgrSolver<D>::printLbmCell(int level, int i, int j, int k, int set) {
+	stdCellId *newcid = new stdCellId;
+	newcid->level = level;
+	newcid->x = i;
+	newcid->y = j;
+	newcid->z = k;
+
+	// this function is not called upon clicking, then its from setMouseClick
+	debugPrintNodeInfo( newcid, set );
+	delete newcid;
+}
+template<class D>
+void 
+LbmFsgrSolver<D>::debugMarkCellCall(int level, int vi,int vj,int vk) {
+	stdCellId *newcid = new stdCellId;
+	newcid->level = level;
+	newcid->x = vi;
+	newcid->y = vj;
+	newcid->z = vk;
+	addCellToMarkedList( newcid );
+}
+
+		
+/*****************************************************************************/
+// implement CellIterator<UniformFsgrCellIdentifier> interface
+/*****************************************************************************/
+
+
+
+// values from guiflkt.cpp
+extern double guiRoiSX, guiRoiSY, guiRoiSZ, guiRoiEX, guiRoiEY, guiRoiEZ;
+extern int guiRoiMaxLev, guiRoiMinLev;
+#define CID_SX (int)( (mLevel[cid->level].lSizex-1) * guiRoiSX )
+#define CID_SY (int)( (mLevel[cid->level].lSizey-1) * guiRoiSY )
+#define CID_SZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiSZ )
+
+#define CID_EX (int)( (mLevel[cid->level].lSizex-1) * guiRoiEX )
+#define CID_EY (int)( (mLevel[cid->level].lSizey-1) * guiRoiEY )
+#define CID_EZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiEZ )
+
+template<class D>
+CellIdentifierInterface* 
+LbmFsgrSolver<D>::getFirstCell( ) {
+	int level = mMaxRefine;
+
+#if LBMDIM==3
+	if(mMaxRefine>0) { level = mMaxRefine-1; } // NO1HIGHESTLEV DEBUG
+#endif
+	level = guiRoiMaxLev;
+	if(level>mMaxRefine) level = mMaxRefine;
+	
+	//errMsg("LbmFsgrSolver::getFirstCell","Celliteration started...");
+	stdCellId *cid = new stdCellId;
+	cid->level = level;
+	cid->x = CID_SX;
+	cid->y = CID_SY;
+	cid->z = CID_SZ;
+	return cid;
+}
+
+template<class D>
+typename LbmFsgrSolver<D>::stdCellId* 
+LbmFsgrSolver<D>::convertBaseCidToStdCid( CellIdentifierInterface* basecid) {
+	//stdCellId *cid = dynamic_cast<stdCellId*>( basecid );
+	stdCellId *cid = (stdCellId*)( basecid );
+	return cid;
+}
+
+template<class D>
+void 
+LbmFsgrSolver<D>::advanceCell( CellIdentifierInterface* basecid) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	if(cid->getEnd()) return;
+
+	//debugOut(" ADb "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10);
+	cid->x++;
+	if(cid->x > CID_EX){ cid->x = CID_SX; cid->y++; 
+		if(cid->y > CID_EY){ cid->y = CID_SY; cid->z++; 
+			if(cid->z > CID_EZ){ 
+				cid->level--;
+				cid->x = CID_SX; 
+				cid->y = CID_SY; 
+				cid->z = CID_SZ; 
+				if(cid->level < guiRoiMinLev) {
+					cid->level = guiRoiMaxLev;
+					cid->setEnd( true );
+				}
+			}
+		}
+	}
+	//debugOut(" ADa "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10);
+}
+
+template<class D>
+bool 
+LbmFsgrSolver<D>::noEndCell( CellIdentifierInterface* basecid) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return (!cid->getEnd());
+}
+
+template<class D>
+void 
+LbmFsgrSolver<D>::deleteCellIterator( CellIdentifierInterface** cid ) {
+	delete *cid;
+	*cid = NULL;
+}
+
+template<class D>
+CellIdentifierInterface* 
+LbmFsgrSolver<D>::getCellAt( ntlVec3Gfx pos ) {
+	//int cellok = false;
+	pos -= (D::mvGeoStart);
+
+	LbmFloat mmaxsize = mLevel[mMaxRefine].nodeSize;
+	for(int level=mMaxRefine; level>=0; level--) { // finest first
+	//for(int level=0; level<=mMaxRefine; level++) { // coarsest first
+		LbmFloat nsize = mLevel[level].nodeSize;
+		int x,y,z;
+		//LbmFloat nsize = getCellSize(NULL)[0]*2.0;
+		x = (int)((pos[0]-0.5*mmaxsize) / nsize );
+		y = (int)((pos[1]-0.5*mmaxsize) / nsize );
+		z = (int)((pos[2]-0.5*mmaxsize) / nsize );
+		if(D::cDimension==2) z = 0;
+
+		// double check...
+		//int level = mMaxRefine;
+		if(x<0) continue;
+		if(y<0) continue;
+		if(z<0) continue;
+		if(x>=mLevel[level].lSizex) continue;
+		if(y>=mLevel[level].lSizey) continue;
+		if(z>=mLevel[level].lSizez) continue;
+
+		// return fluid/if/border cells
+		if( ( (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused)) ) ||
+			  ( (level<mMaxRefine) && (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused|CFEmpty)) ) ) {
+			continue;
+		} // */
+
+		stdCellId *newcid = new stdCellId;
+		newcid->level = level;
+		newcid->x = x;
+		newcid->y = y;
+		newcid->z = z;
+		//errMsg("cellAt",D::mName<<" "<<pos<<" l"<<level<<":"<<x<<","<<y<<","<<z<<" "<<convertCellFlagType2String(RFLAG(level, x,y,z, mLevel[level].setCurr)) );
+		return newcid;
+	}
+
+	return NULL;
+}
+
+
+// INFO functions
+
+template<class D>
+int      
+LbmFsgrSolver<D>::getCellSet      ( CellIdentifierInterface* basecid) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return mLevel[cid->level].setCurr;
+	//return mLevel[cid->level].setOther;
+}
+
+template<class D>
+int      
+LbmFsgrSolver<D>::getCellLevel    ( CellIdentifierInterface* basecid) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return cid->level;
+}
+
+template<class D>
+ntlVec3Gfx   
+LbmFsgrSolver<D>::getCellOrigin   ( CellIdentifierInterface* basecid) {
+	ntlVec3Gfx ret;
+
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	ntlVec3Gfx cs( mLevel[cid->level].nodeSize );
+	if(D::cDimension==2) { cs[2] = 0.0; }
+
+	if(D::cDimension==2) {
+		ret =(D::mvGeoStart -(cs*0.5) + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], (D::mvGeoEnd[2]-D::mvGeoStart[2])*0.5 )
+				+ ntlVec3Gfx(0.0,0.0,cs[1]*-0.25)*cid->level )
+			+getCellSize(basecid);
+	} else {
+		ret =(D::mvGeoStart -(cs*0.5) + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], cid->z *cs[2] ))
+			+getCellSize(basecid);
+	}
+	return (ret);
+}
+
+template<class D>
+ntlVec3Gfx   
+LbmFsgrSolver<D>::getCellSize     ( CellIdentifierInterface* basecid) {
+	// return half size
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	ntlVec3Gfx retvec( mLevel[cid->level].nodeSize * 0.5 );
+	// 2d display as rectangles
+	if(D::cDimension==2) { retvec[2] = 0.0; }
+	return (retvec);
+}
+
+template<class D>
+LbmFloat 
+LbmFsgrSolver<D>::getCellDensity  ( CellIdentifierInterface* basecid,int set) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+
+	LbmFloat rho = 0.0;
+	FORDF0 {
+		rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
+	}
+	return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].stepsize) +1.0; // normal
+	//return ((rho-1.0) * D::mpParam->getCellSize() / D::mpParam->getStepTime()) +1.0;
+}
+
+template<class D>
+LbmVec   
+LbmFsgrSolver<D>::getCellVelocity ( CellIdentifierInterface* basecid,int set) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+
+	LbmFloat ux,uy,uz;
+	ux=uy=uz= 0.0;
+	FORDF0 {
+		ux += D::dfDvecX[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
+		uy += D::dfDvecY[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
+		uz += D::dfDvecZ[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
+	}
+	LbmVec vel(ux,uy,uz);
+	// TODO fix...
+	return (vel * mLevel[cid->level].simCellSize / mLevel[cid->level].stepsize * D::mDebugVelScale); // normal
+}
+
+template<class D>
+LbmFloat   
+LbmFsgrSolver<D>::getCellDf( CellIdentifierInterface* basecid,int set, int dir) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return QCELL(cid->level, cid->x,cid->y,cid->z, set, dir);
+}
+template<class D>
+LbmFloat   
+LbmFsgrSolver<D>::getCellMass( CellIdentifierInterface* basecid,int set) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return QCELL(cid->level, cid->x,cid->y,cid->z, set, dMass);
+}
+template<class D>
+LbmFloat   
+LbmFsgrSolver<D>::getCellFill( CellIdentifierInterface* basecid,int set) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
+	if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFFluid) return 1.0;
+	return 0.0;
+	//return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
+}
+template<class D>
+CellFlagType 
+LbmFsgrSolver<D>::getCellFlag( CellIdentifierInterface* basecid,int set) {
+	stdCellId *cid = convertBaseCidToStdCid(basecid);
+	return RFLAG(cid->level, cid->x,cid->y,cid->z, set);
+}
+
+template<class D>
+LbmFloat 
+LbmFsgrSolver<D>::getEquilDf( int l ) {
+	return D::dfEquil[l];
+}
+
+template<class D>
+int 
+LbmFsgrSolver<D>::getDfNum( ) {
+	return D::cDfNum;
+}
+
+#if LBM_USE_GUI==1
+//! show simulation info (implement SimulationObject pure virtual func)
+template<class D>
+void 
+LbmFsgrSolver<D>::debugDisplay(fluidDispSettings *set){ 
+	//lbmDebugDisplay< LbmFsgrSolver<D> >( set, this ); 
+	lbmDebugDisplay( set ); 
+}
+#endif
+
+/*****************************************************************************/
+// strict debugging functions
+/*****************************************************************************/
+#if FSGR_STRICT_DEBUG==1
+#define STRICT_EXIT *((int *)0)=0;
+
+template<class D>
+int LbmFsgrSolver<D>::debLBMGI(int level, int ii,int ij,int ik, int is) {
+	if(level <  0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
+	if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
+
+	if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	return _LBMGI(level, ii,ij,ik, is);
+};
+
+template<class D>
+CellFlagType& LbmFsgrSolver<D>::debRFLAG(int level, int xx,int yy,int zz,int set){
+	return _RFLAG(level, xx,yy,zz,set);   
+};
+
+template<class D>
+CellFlagType& LbmFsgrSolver<D>::debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) {
+	if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	// warning might access all spatial nbs
+	if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	return _RFLAG_NB(level, xx,yy,zz,set, dir);
+};
+
+template<class D>
+CellFlagType& LbmFsgrSolver<D>::debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) {
+	if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	return _RFLAG_NBINV(level, xx,yy,zz,set, dir);
+};
+
+template<class D>
+int LbmFsgrSolver<D>::debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
+	if(level <  0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
+	if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
+
+	if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
+	if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
+	if(l>D::cDfNum){  // dFfrac is an exception
+		if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
+#if COMPRESSGRIDS==1
+	//if((!D::mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug
+#endif // COMPRESSGRIDS==1
+	return _LBMQI(level, ii,ij,ik, is, l);
+};
+
+template<class D>
+LbmFloat& LbmFsgrSolver<D>::debQCELL(int level, int xx,int yy,int zz,int set,int l) {
+	//errMsg("LbmStrict","debQCELL debug: l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" l"<<l<<" index"<<LBMGI(level, xx,yy,zz,set)); 
+	return _QCELL(level, xx,yy,zz,set,l);
+};
+
+template<class D>
+LbmFloat& LbmFsgrSolver<D>::debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) {
+	if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	return _QCELL_NB(level, xx,yy,zz,set, dir,l);
+};
+
+template<class D>
+LbmFloat& LbmFsgrSolver<D>::debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) {
+	if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
+	return _QCELL_NBINV(level, xx,yy,zz,set, dir,l);
+};
+
+template<class D>
+LbmFloat* LbmFsgrSolver<D>::debRACPNT(int level,  int ii,int ij,int ik, int is ) {
+	return _RACPNT(level, ii,ij,ik, is );
+};
+
+template<class D>
+LbmFloat& LbmFsgrSolver<D>::debRAC(LbmFloat* s,int l) {
+	if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
+	if(l>dTotalNum){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } 
+	//if(l>D::cDfNum){ // dFfrac is an exception 
+	//if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
+	return _RAC(s,l);
+};
+
+#endif // FSGR_STRICT_DEBUG==1
+
+
+
+#if LBM_USE_GUI==1
+#define USE_GLUTILITIES
+#include "../gui/gui_utilities.h"
+
+//! display a single node
+template<class D>
+void LbmFsgrSolver<D>::debugDisplayNode(fluidDispSettings *dispset, CellIdentifierInterface* cell ) {
+	//debugOut(" DD: "<<cell->getAsString() , 10);
+	ntlVec3Gfx org      = this->getCellOrigin( cell );
+	ntlVec3Gfx halfsize = this->getCellSize( cell );
+	int    set      = this->getCellSet( cell );
+	//debugOut(" DD: "<<cell->getAsString()<<" "<< (dispset->type) , 10);
+
+	bool     showcell = true;
+	int      linewidth = 1;
+	ntlColor col(0.5);
+	LbmFloat cscale = dispset->scale;
+
+#define DRAWDISPCUBE(col,scale) \
+	{	glLineWidth( linewidth ); \
+	  glColor3f( (col)[0], (col)[1], (col)[2]); \
+		ntlVec3Gfx s = org-(halfsize * (scale)); \
+		ntlVec3Gfx e = org+(halfsize * (scale)); \
+		drawCubeWire( s,e ); }
+
+	switch(dispset->type) {
+		case FLUIDDISPNothing: {
+				showcell = false;
+			} break;
+		case FLUIDDISPCelltypes: {
+				CellFlagType flag = this->getCellFlag(cell, set );
+				cscale = 0.5;
+
+				if(flag& CFInvalid  ) { if(!guiShowInvalid  ) return; }
+				if(flag& CFUnused   ) { if(!guiShowInvalid  ) return; }
+				if(flag& CFEmpty    ) { if(!guiShowEmpty    ) return; }
+				if(flag& CFInter    ) { if(!guiShowInterface) return; }
+				if(flag& CFNoDelete ) { if(!guiShowNoDelete ) return; }
+				if(flag& CFBnd      ) { if(!guiShowBnd      ) return; }
+
+				// only dismiss one of these types 
+ 				if(flag& CFGrFromCoarse)  { if(!guiShowCoarseInner  ) return; } // inner not really interesting
+				else
+				if(flag& CFGrFromFine) { if(!guiShowCoarseBorder ) return; }
+				else
+				if(flag& CFFluid    )    { if(!guiShowFluid    ) return; }
+
+				if(flag& CFNoDelete) { // debug, mark nodel cells
+					ntlColor ccol(0.7,0.0,0.0);
+					DRAWDISPCUBE(ccol, 0.1);
+				}
+				if(flag& CFPersistMask) { // mark persistent flags
+					ntlColor ccol(0.5);
+					DRAWDISPCUBE(ccol, 0.125);
+				}
+				if(flag& CFNoBndFluid) { // mark persistent flags
+					ntlColor ccol(0,0,1);
+					DRAWDISPCUBE(ccol, 0.075);
+				}
+
+				if(flag& CFInvalid) {
+					cscale = 0.50;
+					col = ntlColor(0.0,0,0.0);
+				}
+				else if(flag& CFBnd) {
+					cscale = 0.59;
+					col = ntlColor(0.4);
+				}
+
+				else if(flag& CFInter) {
+					cscale = 0.55;
+					col = ntlColor(0,1,1);
+
+				} else if(flag& CFGrFromCoarse) {
+					// draw as - with marker
+					ntlColor col2(0.0,1.0,0.3);
+					DRAWDISPCUBE(col2, 0.1);
+					cscale = 0.5;
+					showcell=false; // DEBUG
+				}
+				else if(flag& CFFluid) {
+					cscale = 0.5;
+					if(flag& CFGrToFine) {
+						ntlColor col2(0.5,0.0,0.5);
+						DRAWDISPCUBE(col2, 0.1);
+						col = ntlColor(0,0,1);
+					}
+					if(flag& CFGrFromFine) {
+						ntlColor col2(1.0,1.0,0.0);
+						DRAWDISPCUBE(col2, 0.1);
+						col = ntlColor(0,0,1);
+					} else if(flag& CFGrFromCoarse) {
+						// draw as fluid with marker
+						ntlColor col2(0.0,1.0,0.3);
+						DRAWDISPCUBE(col2, 0.1);
+						col = ntlColor(0,0,1);
+					} else {
+						col = ntlColor(0,0,1);
+					}
+				}
+				else if(flag& CFEmpty) {
+					showcell=false;
+				}
+
+			} break;
+		case FLUIDDISPVelocities: {
+				// dont use cube display
+				LbmVec vel = this->getCellVelocity( cell, set );
+				glBegin(GL_LINES);
+				glColor3f( 0.0,0.0,0.0 );
+				glVertex3f( org[0], org[1], org[2] );
+				org += vec2G(vel * 10.0 * cscale);
+				glColor3f( 1.0,1.0,1.0 );
+				glVertex3f( org[0], org[1], org[2] );
+				glEnd();
+				showcell = false;
+			} break;
+		case FLUIDDISPCellfills: {
+				CellFlagType flag = this->getCellFlag( cell,set );
+				cscale = 0.5;
+
+				if(flag& CFFluid) {
+					cscale = 0.75;
+					col = ntlColor(0,0,0.5);
+				}
+				else if(flag& CFInter) {
+					cscale = 0.75 * this->getCellMass(cell,set);
+					col = ntlColor(0,1,1);
+				}
+				else {
+					showcell=false;
+				}
+
+					if( ABS(this->getCellMass(cell,set)) < 10.0 ) {
+						cscale = 0.75 * this->getCellMass(cell,set);
+					} else {
+						showcell = false;
+					}
+					if(cscale>0.0) {
+						col = ntlColor(0,1,1);
+					} else {
+						col = ntlColor(1,1,0);
+					}
+			// TODO
+			} break;
+		case FLUIDDISPDensity: {
+				LbmFloat rho = this->getCellDensity(cell,set);
+				cscale = rho*rho * 0.25;
+				col = ntlColor( MIN(0.5+cscale,1.0) , MIN(0.0+cscale,1.0), MIN(0.0+cscale,1.0) );
+				cscale *= 2.0;
+			} break;
+		case FLUIDDISPGrid: {
+				cscale = 0.59;
+				col = ntlColor(1.0);
+			} break;
+		default: {
+				cscale = 0.5;
+				col = ntlColor(1.0,0.0,0.0);
+			} break;
+	}
+
+	if(!showcell) return;
+	DRAWDISPCUBE(col, cscale);
+}
+
+//! debug display function
+//  D has to implement the CellIterator interface
+template<class D>
+void LbmFsgrSolver<D>::lbmDebugDisplay(fluidDispSettings *dispset) {
+	//je nach solver...?
+	if(!dispset->on) return;
+	glDisable( GL_LIGHTING ); // dont light lines
+
+	typename D::CellIdentifier cid = this->getFirstCell();
+	for(; this->noEndCell( cid );
+	      this->advanceCell( cid ) ) {
+		this->debugDisplayNode(dispset, cid );
+	}
+	delete cid;
+
+	glEnable( GL_LIGHTING ); // dont light lines
+}
+
+//! debug display function
+//  D has to implement the CellIterator interface
+template<class D>
+void LbmFsgrSolver<D>::lbmMarkedCellDisplay() {
+	fluidDispSettings dispset;
+	// trick - display marked cells as grid displa -> white, big
+	dispset.type = FLUIDDISPGrid;
+	dispset.on = true;
+	glDisable( GL_LIGHTING ); // dont light lines
+	
+	typename D::CellIdentifier cid = this->markedGetFirstCell();
+	while(cid) {
+		this->debugDisplayNode(&dispset, cid );
+		cid = this->markedAdvanceCell();
+	}
+	delete cid;
+
+	glEnable( GL_LIGHTING ); // dont light lines
+}
+
+#endif // LBM_USE_GUI==1
+
+//! display a single node
+template<class D>
+void LbmFsgrSolver<D>::debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet) {
+		//string printInfo,
+		// force printing of one set? default = -1 = off
+  bool printDF     = false;
+  bool printRho    = false;
+  bool printVel    = false;
+  bool printFlag   = false;
+  bool printGeom   = false;
+  bool printMass=false;
+	bool printBothSets = false;
+	string printInfo = this->getNodeInfoString();
+
+	for(size_t i=0; i<printInfo.length()-0; i++) {
+		char what = printInfo[i];
+		switch(what) {
+			case '+': // all on
+								printDF = true; printRho = true; printVel = true; printFlag = true; printGeom = true; printMass = true ;
+								printBothSets = true; break;
+			case '-': // all off
+								printDF = false; printRho = false; printVel = false; printFlag = false; printGeom = false; printMass = false; 
+								printBothSets = false; break;
+			case 'd': printDF = true; break;
+			case 'r': printRho = true; break;
+			case 'v': printVel = true; break;
+			case 'f': printFlag = true; break;
+			case 'g': printGeom = true; break;
+			case 'm': printMass = true; break;
+			case 's': printBothSets = true; break;
+			default: 
+				errFatal("debugPrintNodeInfo","Invalid node info id "<<what,SIMWORLD_GENERICERROR); return;
+		}
+	}
+
+	ntlVec3Gfx org      = this->getCellOrigin( cell );
+	ntlVec3Gfx halfsize = this->getCellSize( cell );
+	int    set      = this->getCellSet( cell );
+	debMsgStd("debugPrintNodeInfo",DM_NOTIFY, "Printing cell info '"<<printInfo<<"' for node: "<<cell->getAsString()<<" from "<<this->getName()<<" currSet:"<<set , 1);
+	if(printGeom) debMsgStd("                  ",DM_MSG, "Org:"<<org<<" Halfsize:"<<halfsize<<" ", 1);
+
+	int setmax = 2;
+	if(!printBothSets) setmax = 1;
+	if(forceSet>=0) setmax = 1;
+
+	for(int s=0; s<setmax; s++) {
+		int workset = set;
+		if(s==1){ workset = (set^1); }		
+		if(forceSet>=0) workset = forceSet;
+		debMsgStd("                  ",DM_MSG, "Printing set:"<<workset<<" orgSet:"<<set, 1);
+		
+		if(printDF) {
+			for(int l=0; l<this->getDfNum(); l++) { // FIXME ??
+				debMsgStd("                  ",DM_MSG, "  Df"<<l<<": "<<this->getCellDf(cell,workset,l), 1);
+			}
+		}
+		if(printRho) {
+			debMsgStd("                  ",DM_MSG, "  Rho: "<<this->getCellDensity(cell,workset), 1);
+		}
+		if(printVel) {
+			debMsgStd("                  ",DM_MSG, "  Vel: "<<this->getCellVelocity(cell,workset), 1);
+		}
+		if(printFlag) {
+			CellFlagType flag = this->getCellFlag(cell,workset);
+			debMsgStd("                  ",DM_MSG, "  Flg: "<< flag<<" "<<convertFlags2String( flag ) <<" "<<convertCellFlagType2String( flag ), 1);
+		}
+		if(printMass) {
+			debMsgStd("                  ",DM_MSG, "  Mss: "<<this->getCellMass(cell,workset), 1);
+		}
+	}
+}
+
+
+
+#if LBMDIM==2
+template class LbmFsgrSolver< LbmBGK2D >;
+#endif // LBMDIM==2
+#if LBMDIM==3
+template class LbmFsgrSolver< LbmBGK3D >;
+#endif // LBMDIM==3

intern/elbeem/intern/typeslbm.h

@@ -1,260 +0,0 @@
-/******************************************************************************
- *
- * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
- * Copyright 2003,2004 Nils Thuerey
- *
- * Lattice-Boltzmann defines...
- * 
- *****************************************************************************/
-#ifndef TYPES_LBM_H
-
-/* standard precision for LBM solver */
-typedef double LBM_Float;
-//typedef float LBM_Float;
-
-typedef double LBM2D_Float;
-//typedef float LBM2D_Float; // FLAGS might not work!!!!
-
-
-/******************************************************************************
- * 2D 
- *****************************************************************************/
-
-
-//! use incompressible LGBK model?
-#define LBM2D_INCOMPBGK 1
-
-
-/*! size of a single set of distribution functions */
-#define LBM2D_DISTFUNCSIZE  9
-/*! size of a single set for a cell (+cell flags, mass, bubble id) */
-#define LBM2D_SETSIZE  12
-/*! floats per LBM cell */
-#define LBM2D_FLOATSPERCELL (LBM2D_SETSIZE +LBM2D_SETSIZE )
-
-
-/*! sphere init full or empty */
-#define LBM2D_FILLED true
-#define LBM2D_EMPTY  false
-
-/*! distribution functions directions */
-#define WC   0
-#define WN   1
-#define WS   2
-#define WE   3
-#define WW   4
-#define WNE  5
-#define WNW  6
-#define WSE  7
-#define WSW  8
-#define FLAG2D_BND    (9)
-#define FLAG2D_MASS   (10)
-#define FLAG2D_BUBBLE (11)
-
-/* Wi factors for collide step */
-#define LBM2D_COLLEN_ZERO   (4.0/9.0)
-#define LBM2D_COLLEN_ONE    (1.0/9.0)
-#define LBM2D_COLLEN_SQRTWO (1.0/36.0)
-
-
-/* calculate equlibrium function for a single direction at cell i,j
- * pass 0 for the u_ terms that are not needed
- */
-#define LBM2D_VELVEC(l, ux,uy) ((ux)*DF2DdvecX[l]+(uy)*DF2DdvecY[l])
-#if LBM2D_INCOMPBGK!=1
-#define LBM2D_COLLIDE_EQ(target, l,Rho, ux,uy)    \
-  {\
-		LBM2D_Float tmp = LBM2D_VELVEC(l,ux,uy); \
-    target = ( (DF2Dlength[l]*Rho) *( \
-		  					+ 1.0 - (3.0/2.0*(ux*ux + uy*uy)) \
-							  + 3.0 *tmp \
-							  + 9.0/2.0 *(tmp*tmp) ) \
-	  );\
-  }
-#endif
-
-/* incompressible LBGK model?? */
-#if LBM2D_INCOMPBGK==1
-#define LBM2D_COLLIDE_EQ(target, l,Rho, ux,uy)    \
-  {\
-		LBM2D_Float tmp = LBM2D_VELVEC(l,ux,uy); \
-    target = ( (DF2Dlength[l]) *( \
-		  					+ Rho - (3.0/2.0*(ux*ux + uy*uy )) \
-							  + 3.0 *tmp \
-							  + 9.0/2.0 *(tmp*tmp) ) \
-	  );\
-  }
-#endif
-
-/* calculate new distribution function for cell i,j
- * Now also includes gravity
- */
-#define LBM2D_COLLIDE(l,omega, Rho, ux,uy ) \
-  {\
-		LBM2D_Float collideTempVar; \
-		LBM2D_COLLIDE_EQ(collideTempVar, l,Rho, (ux), (uy) ); \
-		m[l] = (1.0-omega) * m[l] + \
- 					 omega* collideTempVar \
-           ; \
-  }\
-
-
-#ifdef LBM2D_IMPORT
-extern char *DF2Dstring[LBM2D_DISTFUNCSIZE];
-extern int DF2Dnorm[LBM2D_DISTFUNCSIZE];
-extern int DF2Dinv[LBM2D_DISTFUNCSIZE];
-extern int DF2DrefX[LBM2D_DISTFUNCSIZE];
-extern int DF2DrefY[LBM2D_DISTFUNCSIZE];
-extern LBM2D_Float DF2Dequil[ LBM2D_DISTFUNCSIZE ];
-extern int DF2DvecX[LBM2D_DISTFUNCSIZE];
-extern int DF2DvecY[LBM2D_DISTFUNCSIZE];
-extern LBM2D_Float DF2DdvecX[LBM2D_DISTFUNCSIZE];
-extern LBM2D_Float DF2DdvecY[LBM2D_DISTFUNCSIZE];
-extern LBM2D_Float DF2Dlength[LBM2D_DISTFUNCSIZE];
-#endif
-
-
-
-/******************************************************************************
- * 3D 
- *****************************************************************************/
-
-// use incompressible LGBK model?
-#define LBM_INCOMPBGK 1
-
-
-
-/*! size of a single set of distribution functions */
-#define LBM_DISTFUNCSIZE  19
-/*! size of a single set for a cell (+cell flags, mass, bubble id) */
-#define LBM_SETSIZE  22
-/*! floats per LBM cell */
-#define LBM_FLOATSPERCELL (LBM_SETSIZE +LBM_SETSIZE )
-
-
-/*! distribution functions directions */
-#define MC   0
-#define MN   1
-#define MS   2
-#define ME   3
-#define MW   4
-#define MT   5
-#define MB   6
-#define MNE  7
-#define MNW  8
-#define MSE  9
-#define MSW 10
-#define MNT 11
-#define MNB 12 
-#define MST 13
-#define MSB 14
-#define MET 15 
-#define MEB 16
-#define MWT 17
-#define MWB 18
-#define FLAG_BND    (19)
-#define FLAG_MASS   (20)
-#define FLAG_BUBBLE (21)
-
-/* Wi factors for collide step */
-#define LBM_COLLEN_ZERO   (1.0/3.0)
-#define LBM_COLLEN_ONE    (1.0/18.0)
-#define LBM_COLLEN_SQRTWO (1.0/36.0)
-
-
-/* calculate equlibrium function for a single direction at cell i,j,k
- * pass 0 for the u_ terms that are not needed
- */
-#define LBM_VELVEC(l, ux,uy,uz) ((ux)*DFdvecX[l]+(uy)*DFdvecY[l]+(uz)*DFdvecZ[l])
-#ifndef LBM_INCOMPBGK
-#define LBM_COLLIDE_EQ(target, l,Rho, ux,uy,uz)    \
-  {\
-		LBM_Float tmp = LBM_VELVEC(l,ux,uy,uz); \
-    target = ( (DFlength[l]*Rho) *( \
-		  					+ 1.0 - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) \
-							  + 3.0 *tmp \
-							  + 9.0/2.0 *(tmp*tmp) ) \
-	  );\
-  }
-#endif
-
-/* incompressible LBGK model?? */
-#ifdef LBM_INCOMPBGK
-#define LBM_COLLIDE_EQ(target, l,Rho, ux,uy,uz)    \
-  {\
-		LBM_Float tmp = LBM_VELVEC(l,ux,uy,uz); \
-    target = ( (DFlength[l]) *( \
-		  					+ Rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) \
-							  + 3.0 *tmp \
-							  + 9.0/2.0 *(tmp*tmp) ) \
-	  );\
-  }
-#endif
-
-/* calculate new distribution function for cell i,j,k
- * Now also includes gravity
- */
-#define LBM_COLLIDE(l,omega,Rho, ux,uy,uz ) \
-  {\
-		LBM_Float collideTempVar; \
-		LBM_COLLIDE_EQ(collideTempVar, l,Rho, (ux), (uy), (uz) ); \
-		m[l] = (1.0-omega) * m[l] + \
- 					 omega* collideTempVar \
-           ; \
-  }\
-
-
-#ifdef LBM3D_IMPORT
-char *DFstring[LBM_DISTFUNCSIZE];
-int DFnorm[LBM_DISTFUNCSIZE];
-int DFinv[LBM_DISTFUNCSIZE];
-int DFrefX[LBM_DISTFUNCSIZE];
-int DFrefY[LBM_DISTFUNCSIZE];
-int DFrefZ[LBM_DISTFUNCSIZE];
-LBM_Float DFequil[ LBM_DISTFUNCSIZE ];
-int DFvecX[LBM_DISTFUNCSIZE];
-int DFvecY[LBM_DISTFUNCSIZE];
-int DFvecZ[LBM_DISTFUNCSIZE];
-LBM_Float DFdvecX[LBM_DISTFUNCSIZE];
-LBM_Float DFdvecY[LBM_DISTFUNCSIZE];
-LBM_Float DFdvecZ[LBM_DISTFUNCSIZE];
-LBM_Float DFlength[LBM_DISTFUNCSIZE];
-#endif
-
-
-/******************************************************************************
- * BOTH
- *****************************************************************************/
-
-/*! boundary flags 
- * only 1 should be active for a cell */
-#define BND  (1<< 0)
-#define ACCX (1<< 1)
-#define ACCY (1<< 2)
-#define ACCZ (1<< 3)
-#define FREESLIP (1<< 4)
-#define NOSLIP   (1<< 5)
-#define PERIODIC (1<< 6)
-#define PARTSLIP (1<< 7)
-/*! surface type, also only 1 should be active (2. flag byte) */
-#define EMPTY   (0)
-#define FLUID   (1<< 8)
-#define INTER   (1<< 9)
-/*! neighbor flags (3. flag byte) */
-#define I_NONBFLUID (1<<16)
-#define I_NONBINTER (1<<17)
-#define I_NONBEMPTY (1<<18)
-#define I_NODELETE  (1<<19)
-#define I_NEWCELL   (1<<20)
-#define I_NEWINTERFACE (1<<21)
-/*! marker only for debugging, this bit is reset each step */
-#define I_CELLMARKER     ((int) (1<<30) )
-#define I_NOTCELLMARKER  ((int) (~(1<<30)) )
-
-
-
-
-
-#define TYPES_LBM_H
-#endif
-