- another minor solver update to fix

obstacle fluid surface generation bug
- also contains some code clean ups
  and safer initializations

intern/elbeem/intern/attributes.cpp

@@ -472,35 +472,45 @@ bool AttributeList::ignoreParameter(string name, string source) {
 		
 // read channels
 AnimChannel<int> AttributeList::readChannelInt(string name, int defaultValue, string source, string target, bool needed) {
-	if(!exists(name)) { return AnimChannel<int>(defaultValue); } 
+	if(!exists(name)) { 
+		if(needed) { errFatal("AttributeList::readChannelInt","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+		return AnimChannel<int>(defaultValue); } 
 	AnimChannel<int> ret = find(name)->getChannelInt(); 
 	find(name)->setUsed(true);
 	channelSimplifyi(ret);
 	return ret;
 }
 AnimChannel<double> AttributeList::readChannelFloat(string name, double defaultValue, string source, string target, bool needed ) {
-	if(!exists(name)) { return AnimChannel<double>(defaultValue); } 
+	if(!exists(name)) {
+		if(needed) { errFatal("AttributeList::readChannelFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+	 	return AnimChannel<double>(defaultValue); } 
 	AnimChannel<double> ret = find(name)->getChannelFloat(); 
 	find(name)->setUsed(true);
 	channelSimplifyd(ret);
 	return ret;
 }
 AnimChannel<ntlVec3d> AttributeList::readChannelVec3d(string name, ntlVec3d defaultValue, string source, string target, bool needed ) {
-	if(!exists(name)) { return AnimChannel<ntlVec3d>(defaultValue); } 
+	if(!exists(name)) { 
+		if(needed) { errFatal("AttributeList::readChannelVec3d","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+		return AnimChannel<ntlVec3d>(defaultValue); } 
 	AnimChannel<ntlVec3d> ret = find(name)->getChannelVec3d(); 
 	find(name)->setUsed(true);
 	channelSimplifyVd(ret);
 	return ret;
 }
 AnimChannel<ntlSetVec3f> AttributeList::readChannelSetVec3f(string name, ntlSetVec3f defaultValue, string source, string target, bool needed) {
-	if(!exists(name)) { return AnimChannel<ntlSetVec3f>(defaultValue); } 
+	if(!exists(name)) { 
+		if(needed) { errFatal("AttributeList::readChannelSetVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+		return AnimChannel<ntlSetVec3f>(defaultValue); } 
 	AnimChannel<ntlSetVec3f> ret = find(name)->getChannelSetVec3f(); 
 	find(name)->setUsed(true);
 	//channelSimplifyVf(ret);
 	return ret;
 }
 AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float defaultValue, string source, string target, bool needed ) {
-	if(!exists(name)) { return AnimChannel<float>(defaultValue); } 
+	if(!exists(name)) { 
+		if(needed) { errFatal("AttributeList::readChannelSinglePrecFloat","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+		return AnimChannel<float>(defaultValue); } 
 	AnimChannel<double> convert = find(name)->getChannelFloat(); 
 	find(name)->setUsed(true);
 	channelSimplifyd(convert);
@@ -514,7 +524,9 @@ AnimChannel<float> AttributeList::readChannelSinglePrecFloat(string name, float
 	return ret;
 }
 AnimChannel<ntlVec3f> AttributeList::readChannelVec3f(string name, ntlVec3f defaultValue, string source, string target, bool needed) {
-	if(!exists(name)) { return AnimChannel<ntlVec3f>(defaultValue); } 
+	if(!exists(name)) { 
+		if(needed) { errFatal("AttributeList::readChannelVec3f","Required channel '"<<name<<"' for "<<target<<" from "<< source <<"  not set! ", SIMWORLD_INITERROR); }
+		return AnimChannel<ntlVec3f>(defaultValue); } 
 
 	AnimChannel<ntlVec3d> convert = find(name)->getChannelVec3d(); 
 	// convert to float
@@ -737,6 +749,19 @@ void AnimChannel<Scalar>::debugPrintChannel() {
 }
 
 
+// hack to force instantiation
+void __forceAnimChannelInstantiation() {
+	AnimChannel< float > tmp1;
+	AnimChannel< double > tmp2;
+	AnimChannel< string > tmp3;
+	AnimChannel< ntlVector3Dim<float> > tmp4;
+	tmp1.debugPrintChannel();
+	tmp2.debugPrintChannel();
+	tmp3.debugPrintChannel();
+	tmp4.debugPrintChannel();
+}
+
+
 ntlSetVec3f::ntlSetVec3f(double v ) {
 	mVerts.clear();
 	mVerts.push_back( ntlVec3f(v) );

intern/elbeem/intern/elbeem.cpp

@@ -115,32 +115,45 @@ void elbeemGetErrorString(char *buffer) {
 extern "C" 
 void elbeemResetMesh(elbeemMesh *mesh) {
 	if(!mesh) return;
+	// init typedef struct elbeemMesh
   mesh->type = 0;
-  mesh->parentDomainId = 0;
+
+	mesh->parentDomainId = 0;
+
+	/* vertices */
   mesh->numVertices = 0;
 	mesh->vertices = NULL;
+
+	mesh->channelSizeVertices = 0;
+	mesh->channelVertices = NULL;
+
+	/* triangles */
 	mesh->numTriangles = 0;
   mesh->triangles = NULL;
 
+	/* animation channels */
 	mesh->channelSizeTranslation = 0;
 	mesh->channelTranslation = NULL;
 	mesh->channelSizeRotation = 0;
 	mesh->channelRotation = NULL;
 	mesh->channelSizeScale = 0;
 	mesh->channelScale = NULL;
+
+	/* active channel */
 	mesh->channelSizeActive = 0;
 	mesh->channelActive = NULL;
+
 	mesh->channelSizeInitialVel = 0;
 	mesh->channelInitialVel = NULL;
+
 	mesh->localInivelCoords = 0;
 
 	mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP;
-	mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
 	mesh->obstaclePartslip= 0.;
 
-	mesh->channelSizeVertices = 0;
+	mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
-	mesh->channelVertices = NULL;
 
+	/* name of the mesh, mostly for debugging */
 	mesh->name = "[unnamed]";
 }
 

intern/elbeem/intern/ntl_geometrymodel.cpp

@@ -28,9 +28,7 @@ ntlGeometryObjModel::ntlGeometryObjModel( void ) :
 	mLoaded( false ),
 	mTriangles(), mVertices(), mNormals(),
 	mcAniVerts(), mcAniNorms(), 
-	mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.),
+	mcAniTimes(), mAniTimeScale(1.), mAniTimeOffset(0.)
-	mvCPSStart(-10000.), mvCPSEnd(10000.),
-	mCPSFilename(""), mCPSWidth(0.1), mCPSTimestep(1.)
 {
 }
 
@@ -52,6 +50,39 @@ bool ntlGeometryObjModel::getMeshAnimated() {
 	return ret;
 }
 
+/*! calculate max extends of (ani) mesh */
+void ntlGeometryObjModel::getExtends(ntlVec3Gfx &sstart, ntlVec3Gfx &send) {
+	bool ini=false;
+	ntlVec3Gfx start(0.),end(0.);
+	for(int s=0; s<=(int)mcAniVerts.accessValues().size(); s++) {
+		vector<ntlVec3f> *sverts;
+		if(mcAniVerts.accessValues().size()>0) {
+			if(s==(int)mcAniVerts.accessValues().size()) continue;
+			sverts	= &(mcAniVerts.accessValues()[s].mVerts);
+		} else sverts = &mVertices;
+
+		for(int i=0; i<(int)sverts->size(); i++) {
+
+			if(!ini) {
+				start=(*sverts)[i];
+				end=(*sverts)[i];
+				//errMsg("getExtends","ini "<<s<<","<<i<<" "<<start<<","<<end);
+				ini=true;
+			} else {
+				for(int j=0; j<3; j++) {
+					if(start[j] > (*sverts)[i][j]) { start[j]= (*sverts)[i][j]; }
+					if(end[j]   < (*sverts)[i][j]) { end[j]  = (*sverts)[i][j]; }
+				}
+				//errMsg("getExtends","check "<<s<<","<<i<<" "<<start<<","<<end<<" "<< (*sverts)[i]);
+			}
+
+		}
+	}
+	sstart=start;
+	send=end;
+}
+
+
 /*****************************************************************************/
 /* Init attributes etc. of this object */
 /*****************************************************************************/
@@ -90,11 +121,6 @@ void ntlGeometryObjModel::initialize(ntlRenderGlobals *glob)
 	mAniTimeScale = mpAttrs->readFloat("ani_timescale", mAniTimeScale,"ntlGeometryObjModel", "mAniTimeScale", false);
 	mAniTimeOffset = mpAttrs->readFloat("ani_timeoffset", mAniTimeOffset,"ntlGeometryObjModel", "mAniTimeOffset", false);
 
-	mCPSWidth = mpAttrs->readFloat("cps_width", mCPSWidth,"ntlGeometryObjModel", "mCPSWidth", false);
-	mCPSTimestep = mpAttrs->readFloat("cps_timestep", mCPSTimestep,"ntlGeometryObjModel", "mCPSTimestep", false);
-	mvCPSStart = vec2G( mpAttrs->readVec3d("cps_start", vec2D(mvCPSStart),"ntlGeometryObjModel", "mvCPSStart", false));
-	mvCPSEnd = vec2G( mpAttrs->readVec3d("cps_end", vec2D(mvCPSEnd),"ntlGeometryObjModel", "mvCPSEnd", false));
-
 	// continue with standard obj
 	if(loadBobjModel(mFilename)==0) mLoaded=1;
 	if(!mLoaded) {
@@ -322,7 +348,7 @@ int ntlGeometryObjModel::loadBobjModel(string filename)
 		bytesRead += gzread(gzf, &frameTime, sizeof(frameTime) );
 		//if(bytesRead!=3*sizeof(float)) {
 		if(bytesRead!=sizeof(float)) {
-			debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' no ani sets. ", 10 );
+			debMsgStd("ntlGeometryObjModel::loadBobjModel",DM_MSG, "File '"<<filename<<"' end of gzfile. ", 10 );
 			if(anitimes.size()>0) {
 				// finally init channels and stop reading file
 				mcAniVerts = AnimChannel<ntlSetVec3f>(aniverts,anitimes);
@@ -428,7 +454,7 @@ ntlGeometryObjModel::getTriangles(double t, vector<ntlTriangle> *triangles,
 		(*normals)[startvert+i] = mNormals[i];
 	}
 
-	triangles->reserve(triangles->size() + mTriangles.size() );
+	triangles->reserve(triangles->size() + mTriangles.size()/3 );
 	for(int i=0; i<(int)mTriangles.size(); i+=3) {
 		int trip[3];
 		trip[0] = startvert+mTriangles[i+0];

intern/elbeem/intern/ntl_geometrymodel.h

@@ -44,8 +44,8 @@ class ntlGeometryObjModel : public ntlGeometryObject
 		/*! init triangle divisions */
 		virtual void calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTriangle> &tris, gfxReal fsTri);
 
-		/*! do ani mesh CPS */
+		/*! calculate max extends of (ani) mesh */
-		void calculateCPS(string filename);
+		void getExtends(ntlVec3Gfx &start, ntlVec3Gfx &end);
 
 	private:
 
@@ -71,12 +71,6 @@ class ntlGeometryObjModel : public ntlGeometryObject
 		/*! timing mapping & offset for config files */
 		double mAniTimeScale, mAniTimeOffset;
 
-		/*! ani mesh cps params */
-		ntlVec3Gfx mvCPSStart, mvCPSEnd;
-		string mCPSFilename;
-		gfxReal mCPSWidth, mCPSTimestep;
-
-
 	public:
 
 		/* Access methods */
@@ -87,6 +81,9 @@ class ntlGeometryObjModel : public ntlGeometryObject
 		inline ntlVec3Gfx getEnd( void ){ return mvEnd; }
 		inline void setEnd( const ntlVec3Gfx &set ){ mvEnd = set; }
 
+		inline bool getLoaded( void ){ return mLoaded; }
+		inline void setLoaded( bool set ){ mLoaded = set; }
+
 		/*! set data file name */
 		inline void setFilename(string set) { mFilename = set; }
 };

intern/elbeem/intern/ntl_ray.h

@@ -331,6 +331,8 @@ public:
 		mGeos.push_back( geo ); 
 		geo->setObjectId(mGeos.size());
 	}
+	/*! Add a geo object to the scene, warning - only needed for hand init */
+	inline void addGeoObject(ntlGeometryObject *geo) { mObjects.push_back( geo ); }
 
 	/*! Acces a certain object */
 	inline ntlGeometryObject *getObject(int id) { 

intern/elbeem/intern/particletracer.cpp

@@ -154,7 +154,7 @@ void ParticleTracer::cleanup() {
  *! dump particles if desired 
  *****************************************************************************/
 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) {
-	debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts, 10); // DEBUG
+	debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG
 
 	if(
 			(dumptype==DUMP_FULLGEOMETRY)&&
@@ -305,7 +305,7 @@ void ParticleTracer::checkTrails(double time) {
 /******************************************************************************
  * Get triangles for rendering
  *****************************************************************************/
-void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles, 
+void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, 
 													 vector<ntlVec3Gfx> *vertices, 
 													 vector<ntlVec3Gfx> *normals, int objectId )
 {
@@ -326,7 +326,7 @@ void ParticleTracer::getTriangles(double t, vector<ntlTriangle> *triangles,
 	int segments = mPartSegments;
 	ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2]));
 	ntlVec3Gfx trans = mStart;
-	t = 0.; // doesnt matter
+	time = 0.; // doesnt matter
 
 	for(size_t t=0; t<mPrevs.size()+1; t++) {
 		vector<ParticleObject> *dparts;

intern/elbeem/intern/simulation_object.cpp

@@ -96,6 +96,7 @@ void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) {
 	*mpElbeemSettings = *settings;
 
 	mGeoInitId = settings->domainId+1;
+	debMsgStd("SimulationObject",DM_MSG,"mGeoInitId="<<mGeoInitId<<", domainId="<<settings->domainId, 8);
 }
 
 /******************************************************************************
@@ -115,7 +116,7 @@ int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 	}
 
 
-	this->mGeoInitId = mpAttrs->readInt("geoinitid", this->mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
+	mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
 	//mDimension, mSolverType are deprecated
 	string mSolverType(""); 
 	mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false ); 
@@ -297,6 +298,7 @@ void SimulationObject::step( void )
 		// dont advance for stopped time
 		mpLbm->step();
 		mTime += mpParam->getTimestep();
+//if(mTime>0.001) { errMsg("DEBUG!!!!!!!!","quit mlsu..."); exit(1); } // PROFILE DEBUG TEST!
 	}
 	if(mpLbm->getPanic()) mPanic = true;
 

intern/elbeem/intern/solver_class.h

@@ -106,10 +106,6 @@
 #endif
 #endif
 
-#if LBM_INCLUDE_TESTSOLVERS==1
-#include "solver_test.h"
-#endif // LBM_INCLUDE_TESTSOLVERS==1
-
 /*****************************************************************************/
 /*! cell access classes */
 class UniformFsgrCellIdentifier : 
@@ -467,20 +463,6 @@ class LbmFsgrSolver :
 #		endif // FSGR_STRICT_DEBUG==1
 
 		bool mUseTestdata;
-#if LBM_INCLUDE_TESTSOLVERS==1
-		// test functions
-		LbmTestdata *mpTest;
-		void initTestdata();
-		void destroyTestdata();
-		void handleTestdata();
-		void set3dHeight(int ,int );
-
-		void initCpdata();
-		void handleCpdata();
-	public:
-		// needed for testdata
-		void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz);
-#endif // LBM_INCLUDE_TESTSOLVERS==1
 
 	public: // former LbmModelLBGK  functions
 		// relaxation funtions - implemented together with relax macros

intern/elbeem/intern/solver_init.cpp

@@ -411,7 +411,7 @@ LbmFsgrSolver::~LbmFsgrSolver()
 {
   if(!this->mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; }
 #if COMPRESSGRIDS==1
-	delete mLevel[mMaxRefine].mprsCells[1];
+	delete [] mLevel[mMaxRefine].mprsCells[1];
 	mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
 #endif // COMPRESSGRIDS==1
 
@@ -469,6 +469,12 @@ void LbmFsgrSolver::parseAttrList()
 	// FIXME check needed?
 	mFVArea   = this->mpAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
 
+	// debugging - skip some time...
+	double starttimeskip = 0.;
+	starttimeskip = this->mpAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
+	mSimulationTime += starttimeskip;
+	if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
+
 #if LBM_INCLUDE_TESTSOLVERS==1
 	mUseTestdata = 0;
 	mUseTestdata = this->mpAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
@@ -483,7 +489,7 @@ void LbmFsgrSolver::parseAttrList()
 	mUseTestdata = 0;
 	if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
 #endif // LBM_INCLUDE_TESTSOLVERS!=1
-	if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
+  if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
 
 	if(mMaxRefine==0) mInitialCsmago=0.02;
 	mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
@@ -1184,6 +1190,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
 						if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
 					}
 					break; 
+					// off */
 					/*
 				case FGI_BNDPART: rhomass = BND_FILL;
 					otype = ntype = CFBnd|CFBndPartslip;
@@ -1191,7 +1198,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
 				case FGI_BNDFREE: rhomass = BND_FILL;
 					otype = ntype = CFBnd|CFBndFreeslip;
 					break;
-				// */
+					// off */
 				case FGI_BNDNO:   rhomass = BND_FILL;
 					otype = ntype = CFBnd|CFBndNoslip;
 					break;
@@ -1876,7 +1883,6 @@ void LbmFsgrSolver::initFreeSurfaces() {
 		mLevel[lev].setCurr ^= 1;
 	}
 	// copy back...?
-
 }
 
 /*****************************************************************************/

intern/elbeem/intern/solver_interface.cpp

@@ -44,6 +44,7 @@ LbmSolverInterface::LbmSolverInterface() :
 	mName("lbm_default") ,
 	mpIso( NULL ), mIsoValue(0.499),
 	mSilent(false) , 
+	mLbmInitId(1) ,
 	mpGiTree( NULL ),
 	mpGiObjects( NULL ), mGiObjInside(), mpGlob( NULL ),
 	mRefinementDesired(0),
@@ -284,7 +285,7 @@ void LbmSolverInterface::initGeoTree() {
 	if(mpGiTree != NULL) delete mpGiTree;
 	char treeFlag = (1<<(this->mLbmInitId+4));
 	mpGiTree = new ntlTree( 
-			15, 8,  // warning - fixed values for depth & maxtriangles here...
+			15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
 			scene, treeFlag );
 }
 
@@ -299,6 +300,7 @@ void LbmSolverInterface::freeGeoTree() {
 
 
 int globGeoInitDebug = 0;
+int globGICPIProblems = 0;
 /*****************************************************************************/
 /*! check for a certain flag type at position org */
 bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance) {
@@ -371,7 +373,8 @@ bool LbmSolverInterface::geoInitCheckPointInside(ntlVec3Gfx org, int flags, int
 					if(giObjFirstHistSide[i] !=  1) mess=true;
 				}
 				if(mess) {
-					errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
+					//errMsg("IIIproblem","At "<<org<<" obj "<<i<<" inside:"<<mGiObjInside[i]<<" firstside:"<<giObjFirstHistSide[i] );
+					globGICPIProblems++;
 					mGiObjInside[i]++; // believe first hit side...
 				}
 			}

intern/elbeem/intern/solver_main.cpp

@@ -294,6 +294,8 @@ LbmFsgrSolver::mainLoop(int lev)
 	int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
 	//{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIMcDimension<<" mlsz:"<< mLevel[mMaxRefine].lSizez<<" zmax1:"<<getForZMax1(mMaxRefine) ); } // DEBUG
 #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
+#define LIST_EMPTY(x) mListEmpty.push_back( x );
+#define LIST_FULL(x)  mListFull.push_back( x );
 #endif // PARALLEL==1
 
 	// local to loop
@@ -360,22 +362,15 @@ LbmFsgrSolver::mainLoop(int lev)
 	} // COMPRT
 
 #if PARALLEL==0
-  const int id = 0, Nthrds = 1;
+  //const int id = 0, Nthrds = 1;
-#endif // PARALLEL==1
+	const int jstart = 0;
-  const int Nj = mLevel[mMaxRefine].lSizey;
+	const int jend   = mLevel[mMaxRefine].lSizey;
-	int jstart = 0+( id * (Nj / Nthrds) );
+//#endif // PARALLEL==1
-	int jend   = 0+( (id+1) * (Nj / Nthrds) );
+#else // PARALLEL==1
-  if( ((Nj/Nthrds) *Nthrds) != Nj) {
-    errMsg("LbmFsgrSolver","Invalid domain size Nj="<<Nj<<" Nthrds="<<Nthrds);
-  }
-	// cutoff obstacle boundary
-	if(jstart<1) jstart = 1;
-	if(jend>mLevel[mMaxRefine].lSizey-1) jend = mLevel[mMaxRefine].lSizey-1;
-
-#if PARALLEL==1
 	PARA_INITIALIZE();
-	//errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
+	errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
 #endif // PARALLEL==1
+
   for(int k=kstart;k!=kend;k+=kdir) {
 
 	//errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug
@@ -459,13 +454,14 @@ LbmFsgrSolver::mainLoop(int lev)
 				changeFlag(lev, i,j,k, TSET(lev), CFInter);
 
 				// same as ifemptied for if below
-				LbmPoint emptyp; emptyp.flag = 0;
+				LbmPoint oemptyp; oemptyp.flag = 0;
-				emptyp.x = i; emptyp.y = j; emptyp.z = k;
+				oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
 #if PARALLEL==1
-				calcListEmpty[id].push_back( emptyp );
+				//calcListEmpty[id].push_back( oemptyp );
 #else // PARALLEL==1
-				mListEmpty.push_back( emptyp );
+				//mListEmpty.push_back( oemptyp );
 #endif // PARALLEL==1
+				LIST_EMPTY(oemptyp);
 				calcCellsEmptied++;
 				continue;
 			}
@@ -601,6 +597,7 @@ LbmFsgrSolver::mainLoop(int lev)
 		// for fluid cells - just the f_i difference from streaming to empty cells  ----
 		numRecons = 0;
 		bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
+	//onlyBndnb = false; // DEBUG test off
 
 		FORDF1 { // dfl loop
 			recons[l] = 0;
@@ -1008,6 +1005,7 @@ LbmFsgrSolver::mainLoop(int lev)
 
 		// looks much nicer... LISTTRICK
 #if FSGR_LISTTRICK==1
+		if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
 		if(newFlag&CFNoBndFluid) { // test NEW TEST
 			if(!iffilled) {
 				// remove cells independent from amount of change...
@@ -1035,10 +1033,11 @@ LbmFsgrSolver::mainLoop(int lev)
 			if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
 			filledp.x = i; filledp.y = j; filledp.z = k;
 #if PARALLEL==1
-			calcListFull[id].push_back( filledp );
+			//calcListFull[id].push_back( filledp );
 #else // PARALLEL==1
-			mListFull.push_back( filledp );
+			//mListFull.push_back( filledp );
 #endif // PARALLEL==1
+			LIST_FULL(filledp);
 			//this->mNumFilledCells++; // DEBUG
 			calcCellsFilled++;
 		}
@@ -1047,10 +1046,11 @@ LbmFsgrSolver::mainLoop(int lev)
 			if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
 			emptyp.x = i; emptyp.y = j; emptyp.z = k;
 #if PARALLEL==1
-			calcListEmpty[id].push_back( emptyp );
+			//calcListEmpty[id].push_back( emptyp );
 #else // PARALLEL==1
-			mListEmpty.push_back( emptyp );
+			//mListEmpty.push_back( emptyp );
 #endif // PARALLEL==1
+			LIST_EMPTY(emptyp);
 			//this->mNumEmptiedCells++; // DEBUG
 			calcCellsEmptied++;
 		} 

intern/elbeem/intern/solver_relax.h

@@ -17,7 +17,6 @@
 
 
 
-#if LBM_INCLUDE_TESTSOLVERS!=1
 
 // off for non testing
 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
@@ -25,11 +24,8 @@
 	uy += (grav)[1]; \
 	uz += (grav)[2]; \
 
-#else // LBM_INCLUDE_TESTSOLVERS!=1
+#define TEST_IF_CHECK 
 
-// defined in test.h
-
-#endif // LBM_INCLUDE_TESTSOLVERS!=1
 
 	
 /******************************************************************************
@@ -1125,7 +1121,7 @@ inline void LbmFsgrSolver::collideArrays(
 	for(l=0; l<this->cDfNum; l++) { 
 		df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
 	}  
-	if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df);
+	//if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df);
 
 	mux = ux;
 	muy = uy;

intern/elbeem/intern/solver_util.cpp

@@ -53,33 +53,60 @@ void LbmFsgrSolver::prepareVisualization( void ) {
    for(int j=1;j<mLevel[lev].lSizey-1;j++) 
     for(int i=1;i<mLevel[lev].lSizex-1;i++) {
 			const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
-			//continue; // OFF DEBUG
+			//if(cflag&(CFBnd|CFEmpty)) {
-			if(cflag&(CFBnd|CFEmpty)) {
+			if(cflag&(CFBnd)) {
 				continue;
 
-			} else if( (cflag&CFInter) ) {
+			} else if( (cflag&CFEmpty) ) {
-				//} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) && (cflag&CFNoNbFluid) ) {
+				//continue; // OFF DEBUG
-				//} else if( (cflag&CFInter) && (!(cflag&CFNoBndFluid)) ) {
+				int noslipbnd = 0;
+				int intercnt = 0;
+				CellFlagType nbored;
+				FORDF1 { 
+					const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
+					if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
+					if(nbflag&CFInter){ intercnt++; }
+					nbored |= nbflag;
+				}
+				//? val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
+				if((noslipbnd)&&(intercnt>6)) {
+					//if(val<minval) val = minval; 
+					//*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
+					*this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
+				} else if((noslipbnd)&&(intercnt>0)) {
+					*this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.95;
+				} else {
+				}
+				continue;
+
+			//} else if( (cflag&CFInter) ) {
+
+			} else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
+				// optimized fluid
+				val = 1.;
+
+			} else if( (cflag&(CFInter|CFFluid)) ) {
 				int noslipbnd = 0;
 				FORDF1 { 
 					const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
 					if((nbflag&CFBnd)&&(nbflag&CFBnd)&&(CFBndNoslip)){ noslipbnd=1; l=100; } 
 				}
+				// no empty nb interface cells are treated as full
+				if(cflag&(CFNoNbEmpty|CFFluid)) {
+					val=1.0;
+				}
 				val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
+				
 				if(noslipbnd) {
 					//errMsg("NEWVAL", PRINT_IJK<<" val"<<val <<" f"<< convertCellFlagType2String(cflag)<<" "<<noslipbnd); //" nbfl"<<convertCellFlagType2String(nbored) );
 					if(val<minval) val = minval; 
-					*this->mpIso->lbmGetData( i   , j    ,ZKOFF  ) += minval-( val * mIsoWeight[13] ); 
+					*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
 				}
 				// */
-				
-				// no empty nb interface cells are treated as full
-				if(cflag&CFNoNbEmpty) {
-					val=1.0;
-				}
 
-			} else { // fluid?
+			} else { // unused?
-				val = 1.0; 
+				continue;
+				// old fluid val = 1.0; 
 			} // */
 
 			*this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
@@ -885,10 +912,11 @@ void LbmFsgrSolver::advanceParticles() {
 }
 
 void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
+	int workSet = mLevel[mMaxRefine].setCurr;
+	std::ostringstream name;
+
 	// debug - raw dump of ffrac values
 	if(getenv("ELBEEM_RAWDEBUGDUMP")) {
-		int workSet = mLevel[mMaxRefine].setCurr;
-		std::ostringstream name;
 		//name <<"fill_" << this->mStepCnt <<".dump";
 		name << outfilename<< frameNrStr <<".dump";
 		FILE *file = fopen(name.str().c_str(),"w");
@@ -898,9 +926,12 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt
 				for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
 					for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
 						float val = 0.;
-						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
+							val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+							if(val<0.) val=0.;
+							if(val>1.) val=1.;
+						}
 						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
-						//fwrite( &val, sizeof(val), 1, file); // binary
 						fprintf(file, "%f ",val); // text
 						//errMsg("W", PRINT_IJK<<" val:"<<val);
 					}
@@ -912,6 +943,27 @@ void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrSt
 
 		} // file
 	} // */
+	if(getenv("ELBEEM_BINDEBUGDUMP")) {
+		name << outfilename<< frameNrStr <<".bdump";
+		FILE *file = fopen(name.str().c_str(),"w");
+		if(file) {
+			for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
+				for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
+					for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
+						float val = 0.;
+						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
+							val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
+							if(val<0.) val=0.;
+							if(val>1.) val=1.;
+						}
+						if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
+						fwrite( &val, sizeof(val), 1, file); // binary
+					}
+				}
+			}
+			fclose(file);
+		} // file
+	}
 
 	dumptype = 0; frameNr = 0; // get rid of warning
 }