elbeem update:
authorNils Thuerey <nils@thuerey.de>
Tue, 22 Aug 2006 10:45:19 +0000 (10:45 +0000)
committerNils Thuerey <nils@thuerey.de>
Tue, 22 Aug 2006 10:45:19 +0000 (10:45 +0000)
- slightly improved and optimized
  handling of moving obstacles
- cleanup of debug messages and minor fixes

15 files changed:
intern/elbeem/intern/elbeem.h
intern/elbeem/intern/isosurface.cpp
intern/elbeem/intern/isosurface.h
intern/elbeem/intern/loop_tools.h [new file with mode: 0644]
intern/elbeem/intern/ntl_geometryobject.cpp
intern/elbeem/intern/ntl_geometryobject.h
intern/elbeem/intern/simulation_object.cpp
intern/elbeem/intern/solver_adap.cpp
intern/elbeem/intern/solver_class.h
intern/elbeem/intern/solver_init.cpp
intern/elbeem/intern/solver_interface.cpp
intern/elbeem/intern/solver_interface.h
intern/elbeem/intern/solver_main.cpp
intern/elbeem/intern/solver_relax.h
intern/elbeem/intern/solver_util.cpp

index 8eb8a5433b195cc6643a97f14f0c267a6d44cc5a..ea21e807a6f9b45936642eaa5bc9bc33a99fcfef 100644 (file)
@@ -78,6 +78,7 @@ typedef struct elbeemSimulationSettings {
        short generateVertexVectors;
        /* strength of surface smoothing */
        float surfaceSmoothing;
+       // TODO add surf gen flags
 
        /* global transformation to apply to fluidsim mesh */
        float surfaceTrafo[4*4];
index 373a6f747610c2c9d08a60639d543b48fc9aaa32..283fd17eb944122a6391fc53481a5cc0a18d25ae 100644 (file)
@@ -86,6 +86,12 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
 
 
 
+/*! Reset all values */
+void IsoSurface::resetAll(gfxReal val) {
+       int nodes = mSizez*mSizey*mSizex;
+  for(int i=0;i<nodes;i++) { mpData[i] = val; }
+}
+
 
 /******************************************************************************
  * Destructor
@@ -174,6 +180,10 @@ void IsoSurface::triangulate( void )
                                value[6] = *getData(i+1,j+1,k+1);
                                value[7] = *getData(i  ,j+1,k+1);
 
+                               /*int bndskip = 0; // BNDOFFT
+                               for(int s=0; s<8; s++) if(value[s]==-76.) bndskip++;
+                               if(bndskip>0) continue; // */
+
                                // check intersections of isosurface with edges, and calculate cubie index
                                cubeIndex = 0;
                                if (value[0] < mIsoValue) cubeIndex |= 1;
index d1ff15290693b024dcbb4e21269328a4fb005958..5b5198171f354d29d3ecc3e6b68092441f1f64bf 100644 (file)
@@ -47,6 +47,9 @@ class IsoSurface :
                /*! Init ararys etc. */
                virtual void initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent);
 
+               /*! Reset all values */
+               void resetAll(gfxReal val);
+
                /*! triangulate the scalar field given by pointer*/
                void triangulate( void );
 
diff --git a/intern/elbeem/intern/loop_tools.h b/intern/elbeem/intern/loop_tools.h
new file mode 100644 (file)
index 0000000..927263b
--- /dev/null
@@ -0,0 +1,105 @@
+
+// advance pointer in main loop
+#define ADVANCE_POINTERS(p)    \
+       ccel += (QCELLSTEP*(p));        \
+       tcel += (QCELLSTEP*(p));        \
+       pFlagSrc+= (p); \
+       pFlagDst+= (p); \
+       i+= (p);
+
+#define MAX_CALC_ARR 4
+
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+// init region vars
+#define  GRID_REGION_INIT()   \
+       const int istart = -1+gridLoopBound; \
+       const int iend   = mLevel[mMaxRefine].lSizex-1-gridLoopBound; \
+       LbmFloat calcCurrentMass=0; \
+       LbmFloat calcCurrentVolume=0; \
+       int      calcCellsFilled=0; \
+       int      calcCellsEmptied=0; \
+       int      calcNumUsedCells=0; \
+
+
+
+
+//  -----------------------------------------------------------------------------------
+// serial stuff
+#if PARALLEL!=1
+
+#define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
+#define LIST_EMPTY(x) mListEmpty.push_back( x );
+#define LIST_FULL(x)  mListFull.push_back( x );
+
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+#define  GRID_REGION_START()  \
+       { /* main_region */ \
+       int kstart=getForZMinBnd(), kend=getForZMaxBnd(mMaxRefine); \
+       if(gridLoopBound>0){ kstart=getForZMin1(), kend=getForZMax1(mMaxRefine); } \
+       int kdir = 1; \
+       int jstart = gridLoopBound; \
+       int jend   = mLevel[mMaxRefine].lSizey-gridLoopBound; \
+       const int id=0; \
+       LbmFloat *ccel = NULL, *tcel = NULL; \
+       CellFlagType *pFlagSrc=NULL, *pFlagDst=NULL; \
+       if(mLevel[mMaxRefine].setCurr==1) { \
+       kdir = -1; \
+       int temp = kend; \
+       kend = kstart-1; \
+       kstart = temp-1; \
+       temp = id; /* dummy remove warning */ \
+       } \
+
+
+
+       
+#define unused_GRID_REGION_END() \
+       } /* main_region */  \
+       // end unusedGRID_REGION_END
+
+
+//  -----------------------------------------------------------------------------------
+#else // PARALLEL==1
+
+#include "paraloop.h"
+
+#endif // PARALLEL==1
+
+
+//  -----------------------------------------------------------------------------------
+
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+#define  GRID_LOOP_START()   \
+       for(int k=kstart;k!=kend;k+=kdir) { \
+       pFlagSrc = &RFLAG(lev, istart, jstart, k, SRCS(lev)); \
+       pFlagDst = &RFLAG(lev, istart, jstart, k, TSET(lev)); \
+       ccel = RACPNT(lev,     istart, jstart, k, SRCS(lev)); \
+       tcel = RACPNT(lev,     istart, jstart, k, TSET(lev)); \
+       for(int j=jstart;j!=jend;++j) { \
+       /* for(int i=0;i<mLevel[lev].lSizex-2;   ) { */ \
+       for(int i=istart;i!=iend;   ) { \
+       ADVANCE_POINTERS(1); \
+
+
+
+
+// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+#define  GRID_LOOPREG_END()  \
+        \
+       } /* i */ \
+       int i=0; \
+       ADVANCE_POINTERS(2*gridLoopBound); \
+       } /* j */ \
+       } /* all cell loop k,j,i */ \
+       if(doReduce) { } /* dummy remove warning */ \
+       } /* main_region */ \
+        \
+
+
+
+// old loop for COMPRESSGRIDS==0
+#define old__GRID_LOOP_START() \
+  for(int k=kstart;k<kend;++k) { \
+         for(int j=1;j<mLevel[lev].lSizey-1;++j) { \
+               for(int i=0;i<mLevel[lev].lSizex-2;   ) {
+
index 90c3b3437eadfa47e00f464aa020695fa2bb4cea..d815fb55287936d1c08a06d67bf20e39f02c0a2e 100644 (file)
@@ -33,9 +33,9 @@ ntlGeometryObject::ntlGeometryObject() :
        mInitialPos(0.),
        mcTrans(0.), mcRot(0.), mcScale(1.),
        mIsAnimated(false),
-       mMovPoints(), //mMovNormals(),
+       mMovPoints(), mMovNormals(),
        mHaveCachedMov(false),
-       mCachedMovPoints(), //mCachedMovNormals(),
+       mCachedMovPoints(), mCachedMovNormals(),
        mTriangleDivs1(), mTriangleDivs2(), mTriangleDivs3(),
        mMovPntsInited(-100.0), mMaxMovPnt(-1),
        mcGeoActive(1.)
@@ -169,14 +169,6 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
        // always use channel
        if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel<double>(geoactive); }
 
-       //if(    (mcTrans.accessValues().size()>1)  // VALIDATE
-           //|| (mcRot.accessValues().size()>1) 
-           //|| (mcScale.accessValues().size()>1) 
-           //|| (mcGeoActive.accessValues().size()>1) 
-           //|| (mcInitialVelocity.accessValues().size()>1) 
-               //) {
-               //mIsAnimated = true;
-       //}
        checkIsAnimated();
 
        mIsInitialized = true;
@@ -340,14 +332,6 @@ void ntlGeometryObject::initChannels(
        if((act)&&(nAct>0)) {      ADD_CHANNEL_FLOAT(mcGeoActive, nAct, act); }
        if((ivel)&&(nIvel>0)) {    ADD_CHANNEL_VEC(mcInitialVelocity, nIvel, ivel); }
 
-       //if(    (mcTrans.accessValues().size()>1)  // VALIDATE
-           //|| (mcRot.accessValues().size()>1) 
-           //|| (mcScale.accessValues().size()>1) 
-           //|| (mcGeoActive.accessValues().size()>1) 
-           //|| (mcInitialVelocity.accessValues().size()>1) 
-               //) {
-               //mIsAnimated = true;
-       //}
        checkIsAnimated();
        if(debugInitc) { 
                debMsgStd("ntlGeometryObject::initChannels",DM_MSG,getName()<<
@@ -422,6 +406,9 @@ void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTr
        mTriangleDivs1.resize( tris.size() );
        mTriangleDivs2.resize( tris.size() );
        mTriangleDivs3.resize( tris.size() );
+
+       //fsTri *= 2.; // DEBUG! , wrong init!
+
        for(size_t i=0; i<tris.size(); i++) {
                const ntlVec3Gfx p0 = verts[ tris[i].getPoints()[0] ];
                const ntlVec3Gfx p1 = verts[ tris[i].getPoints()[1] ];
@@ -432,17 +419,7 @@ void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTr
                int divs1=0, divs2=0, divs3=0;
                if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
                if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
-               //if(normNoSqrt(side3) > fsTri*fsTri) { divs3 = (int)(norm(side3)/fsTri); }
-               /*if(getMeshAnimated()) {
-                       vector<ntlSetVec3f> *verts =mcAniVerts.accessValues();
-                       for(int s=0; s<verts->size(); s++) {
-                               int tdivs1=0, tdivs2=0, tdivs3=0;
-                               if(normNoSqrt(side1) > fsTri*fsTri) { tdivs1 = (int)(norm(side1)/fsTri); }
-                               if(normNoSqrt(side2) > fsTri*fsTri) { tdivs2 = (int)(norm(side2)/fsTri); }
-                               if(tdivs1>divs1) divs1=tdivs1;
-                               if(tdivs2>divs2) divs2=tdivs2;
-                       }
-               }*/
+
                mTriangleDivs1[i] = divs1;
                mTriangleDivs2[i] = divs2;
                mTriangleDivs3[i] = divs3;
@@ -454,15 +431,14 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
        if((mMovPntsInited==featureSize)&&(!getMeshAnimated())) return;
        const bool debugMoinit=false;
 
-       //vector<ntlVec3Gfx> movNormals;
        vector<ntlTriangle> triangles; 
        vector<ntlVec3Gfx> vertices; 
        vector<ntlVec3Gfx> vnormals; 
        int objectId = 1;
        this->getTriangles(time, &triangles,&vertices,&vnormals,objectId);
        
-       mMovPoints.clear(); //= vertices;
-       //movNormals.clear(); //= vnormals;
+       mMovPoints.clear();
+       mMovNormals.clear();
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" has v:"<<vertices.size()<<" t:"<<triangles.size() );
        // no points?
        if(vertices.size()<1) {
@@ -523,7 +499,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
                        if(dot(mInitialVelocity,n)<0.0) continue;
                }
                mMovPoints.push_back(p);
-               //movNormals.push_back(n);
+               mMovNormals.push_back(n);
                //errMsg("ntlGeometryObject::initMovingPoints","std"<<i<<" p"<<p<<" n"<<n<<" ");
        }
        // init points & refine...
@@ -533,9 +509,9 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
                const ntlVec3Gfx side1 = vertices[ trips[1] ] - p0;
                const ntlVec3Gfx side2 = vertices[ trips[2] ] - p0;
                int divs1=mTriangleDivs1[i], divs2=mTriangleDivs2[i];
-               //if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
-               //if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
-               const ntlVec3Gfx trinorm = getNormalized(cross(side1,side2))*0.5*featureSize;
+               
+               const ntlVec3Gfx trinormOrg = getNormalized(cross(side1,side2));
+               const ntlVec3Gfx trinorm = trinormOrg*0.25*featureSize;
                if(discardInflowBack) { 
                        if(dot(mInitialVelocity,trinorm)<0.0) continue;
                }
@@ -557,25 +533,16 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
                                        //normalize(n);
                                        // discard inflow backsides
 
-                                       mMovPoints.push_back(p);
+                                       mMovPoints.push_back(p + trinorm); // NEW!?
                                        mMovPoints.push_back(p - trinorm);
-                                       //movNormals.push_back(n);
-                                       //movNormals.push_back(trinorm);
+                                       mMovNormals.push_back(trinormOrg);
+                                       mMovNormals.push_back(trinormOrg); 
                                        //errMsg("TRINORM","p"<<p<<" n"<<n<<" trin"<<trinorm);
                                }
                        }
                }
        }
 
-       // duplicate insides
-       //size_t mpsize = mMovPoints.size();
-       //for(size_t i=0; i<mpsize; i++) {
-               //normalize(vnormals[i]);
-               //errMsg("TTAT"," moved:"<<(mMovPoints[i] - mMovPoints[i]*featureSize)<<" org"<<mMovPoints[i]<<" norm"<<mMovPoints[i]<<" fs"<<featureSize);
-               //mMovPoints.push_back(mMovPoints[i] - movNormals[i]*0.5*featureSize);
-               //movNormals.push_back(movNormals[i]);
-       //}
-
        // find max point
        mMaxMovPnt = 0;
        gfxReal dist = normNoSqrt(mMovPoints[0]);
@@ -594,7 +561,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
                // also do trafo...
        } else {
                mCachedMovPoints = mMovPoints;
-               //mCachedMovNormals = movNormals;
+               mCachedMovNormals = mMovNormals;
                //applyTransformation(time, &mCachedMovPoints, &mCachedMovNormals, 0, mCachedMovPoints.size(), true);
                applyTransformation(time, &mCachedMovPoints, NULL, 0, mCachedMovPoints.size(), true);
                mHaveCachedMov = true;
@@ -610,31 +577,27 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
 void ntlGeometryObject::initMovingPointsAnim(
                 double srctime, vector<ntlVec3Gfx> &srcmovPoints,
                 double dsttime, vector<ntlVec3Gfx> &dstmovPoints,
+                vector<ntlVec3Gfx> *dstmovNormals,
                 gfxReal featureSize,
                 ntlVec3Gfx geostart, ntlVec3Gfx geoend
                 ) {
        const bool debugMoinit=false;
 
-       //vector<ntlVec3Gfx> srcmovNormals;
-       //vector<ntlVec3Gfx> dstmovNormals;
-
        vector<ntlTriangle> srctriangles; 
        vector<ntlVec3Gfx> srcvertices; 
        vector<ntlVec3Gfx> unused_normals; 
        vector<ntlTriangle> dsttriangles; 
        vector<ntlVec3Gfx> dstvertices; 
-       //vector<ntlVec3Gfx> dstnormals; 
+       vector<ntlVec3Gfx> dstnormals; 
        int objectId = 1;
        // TODO optimize? , get rid of normals?
        unused_normals.clear();
        this->getTriangles(srctime, &srctriangles,&srcvertices,&unused_normals,objectId);
        unused_normals.clear();
-       this->getTriangles(dsttime, &dsttriangles,&dstvertices,&unused_normals,objectId);
+       this->getTriangles(dsttime, &dsttriangles,&dstvertices,&dstnormals,objectId);
        
        srcmovPoints.clear();
        dstmovPoints.clear();
-       //srcmovNormals.clear();
-       //dstmovNormals.clear();
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has srcv:"<<srcvertices.size()<<" srct:"<<srctriangles.size() );
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has dstv:"<<dstvertices.size()<<" dstt:"<<dsttriangles.size() );
        // no points?
@@ -668,7 +631,7 @@ void ntlGeometryObject::initMovingPointsAnim(
        }
        for(size_t i=0; i<dstvertices.size(); i++) {
                dstmovPoints.push_back(dstvertices[i]);
-               //dstmovNormals.push_back(dstnormals[i]);
+               if(dstmovNormals) (*dstmovNormals).push_back(dstnormals[i]);
        }
        if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<" dst:"<<dstmovPoints.size()<<" " );
        // init points & refine...
@@ -684,8 +647,9 @@ void ntlGeometryObject::initMovingPointsAnim(
                        const ntlVec3Gfx dstp0 =    dstvertices[ dsttrips[0] ];
                        const ntlVec3Gfx dstside1 = dstvertices[ dsttrips[1] ] - dstp0;
                        const ntlVec3Gfx dstside2 = dstvertices[ dsttrips[2] ] - dstp0;
-                       const ntlVec3Gfx src_trinorm = getNormalized(cross(srcside1,srcside2))*0.5*featureSize;
-                       const ntlVec3Gfx dst_trinorm = getNormalized(cross(dstside1,dstside2))*0.5*featureSize;
+                       const ntlVec3Gfx src_trinorm    = getNormalized(cross(srcside1,srcside2))*0.25*featureSize;
+                       const ntlVec3Gfx dst_trinormOrg = getNormalized(cross(dstside1,dstside2));
+                       const ntlVec3Gfx dst_trinorm    = dst_trinormOrg                         *0.25*featureSize;
                        //errMsg("ntlGeometryObject::initMovingPointsAnim","Tri1 "<<srcvertices[srctrips[0]]<<","<<srcvertices[srctrips[1]]<<","<<srcvertices[srctrips[2]]<<" "<<divs1<<","<<divs2 );
                        for(int u=0; u<=divs1; u++) {
                                for(int v=0; v<=divs2; v++) {
@@ -709,60 +673,19 @@ void ntlGeometryObject::initMovingPointsAnim(
                                        if((srcp[1]>geoend[1]  ) && (dstp[1]>geoend[1]  )) continue;
                                        if((srcp[2]>geoend[2]  ) && (dstp[2]>geoend[2]  )) continue;
                                        
-                                       //ntlVec3Gfx srcn = 
-                                               //srcnormals[ srctriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
-                                               //srcnormals[ srctriangles[i].getPoints()[1] ] * uf +
-                                               //srcnormals[ srctriangles[i].getPoints()[2] ] * vf;
-                                       //normalize(srcn);
-                                       srcmovPoints.push_back(srcp);
+                                       srcmovPoints.push_back(srcp+src_trinorm); // SURFENHTEST
                                        srcmovPoints.push_back(srcp-src_trinorm);
-                                       //srcmovNormals.push_back(srcn);
-
-                                       //ntlVec3Gfx dstn = 
-                                               //dstnormals[ dsttriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
-                                               //dstnormals[ dsttriangles[i].getPoints()[1] ] * uf +
-                                               //dstnormals[ dsttriangles[i].getPoints()[2] ] * vf;
-                                       //normalize(dstn);
-                                       dstmovPoints.push_back(dstp);
+
+                                       dstmovPoints.push_back(dstp+dst_trinorm); // SURFENHTEST
                                        dstmovPoints.push_back(dstp-dst_trinorm);
-                                       //dstmovNormals.push_back(dstn);
-
-                                       /*if(debugMoinit && (i>=0)) errMsg("ntlGeometryObject::initMovingPointsAnim"," "<<
-                                       //srcmovPoints[ srcmovPoints.size()-1 ]<<","<<
-                                       //srcmovNormals[ srcmovNormals.size()-1 ]<<"   "<<
-                                       //dstmovPoints[ dstmovPoints.size()-1 ]<<","<<
-                                       //dstmovNormals[ dstmovNormals.size()-1 ]<<" " 
-                                       (srcmovNormals[ srcmovPoints.size()-1 ]-
-                                       dstmovNormals[ dstmovPoints.size()-1 ])<<" "
-                                       );
-                                       // */
+                                       if(dstmovNormals) {
+                                               (*dstmovNormals).push_back(dst_trinormOrg);
+                                               (*dstmovNormals).push_back(dst_trinormOrg); }
                                }
                        }
                }
        }
 
-       /*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
-       // duplicate insides
-       size_t mpsize = srcmovPoints.size();
-       for(size_t i=0; i<mpsize; i++) {
-               srcmovPoints.push_back(srcmovPoints[i] - srcmovNormals[i]*1.0*featureSize);
-               //? srcnormals.push_back(srcnormals[i]);
-       }
-       mpsize = dstmovPoints.size();
-       for(size_t i=0; i<mpsize; i++) {
-               dstmovPoints.push_back(dstmovPoints[i] - dstmovNormals[i]*1.0*featureSize);
-               //? dstnormals.push_back(dstnormals[i]);
-       }
-       // */
-
-       /*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
-       for(size_t i=0; i<srcmovPoints.size(); i++) {
-               ntlVec3Gfx p1 = srcmovPoints[i];
-               ntlVec3Gfx p2 = dstmovPoints[i];
-         gfxReal len = norm(p1-p2);
-               if(len>0.01) { errMsg("ntlGeometryObject::initMovingPointsAnim"," i"<<i<<" "<< p1<<" "<<p2<<" "<<len); }
-       } // */
-
        // find max point not necessary
        debMsgStd("ntlGeometryObject::initMovingPointsAnim",DM_MSG,"Object "<<getName()<<" inited v:"<<srcvertices.size()<<"->"<<srcmovPoints.size()<<","<<dstmovPoints.size() , 5);
 }
@@ -772,8 +695,8 @@ void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3G
        if(mHaveCachedMov) {
                ret = mCachedMovPoints;
                if(norms) { 
-                       //*norms = mCachedMovNormals; 
-                       errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
+                       *norms = mCachedMovNormals; 
+                       //errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
                }
                //errMsg ("ntlGeometryObject::getMovingPoints","Object "<<getName()<<" used cached points "); // DEBUG
                return;
@@ -781,8 +704,8 @@ void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3G
 
        ret = mMovPoints;
        if(norms) { 
-               errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
-               //*norms = mMovNormals; 
+               //errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
+               *norms = mMovNormals; 
        }
 }
 
index 5fe994682e3d5456cf89fb4291447db40766be6c..a9a8109ba1b703930d4b716f3d92231d1a08afa8 100644 (file)
@@ -118,6 +118,7 @@ class ntlGeometryObject : public ntlGeometryClass
                void initMovingPointsAnim(
                 double srctime, vector<ntlVec3Gfx> &srcpoints,
                 double dsttime, vector<ntlVec3Gfx> &dstpoints,
+                vector<ntlVec3Gfx> *dstnormals,
                 gfxReal featureSize, ntlVec3Gfx geostart, ntlVec3Gfx geoend );
                /*! Prepare points for moving objects (copy into ret) */
                void getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3Gfx> *norms = NULL);
@@ -181,11 +182,11 @@ class ntlGeometryObject : public ntlGeometryClass
                
                /*! moving point/normal storage */
                vector<ntlVec3Gfx> mMovPoints;
-               // unused vector<ntlVec3Gfx> mMovNormals;
+               vector<ntlVec3Gfx> mMovNormals;
                /*! cached points for non moving objects/timeslots */
                bool mHaveCachedMov;
                vector<ntlVec3Gfx> mCachedMovPoints;
-               // unused vector<ntlVec3Gfx> mCachedMovNormals;
+               vector<ntlVec3Gfx> mCachedMovNormals;
                /*! precomputed triangle divisions */
                vector<int> mTriangleDivs1,mTriangleDivs2,mTriangleDivs3;
                /*! inited? */
index e34030180a59300dc3f2592b3b5f887c0071973b..3f954f1a4c77ccc1bf0205f9c241419972b55550 100644 (file)
@@ -439,7 +439,7 @@ int SimulationObject::checkCallerStatus(int status, int frame) {
                }
        }
 
-       errMsg("SimulationObject::checkCallerStatus","s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
+       //debMsgStd("SimulationObject::checkCallerStatus",DM_MSG, "s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
        if(isSimworldOk()) return 0;
        return 1;
 }
index aade9b30545cb63be584b086e92227321494cf51..99892081fca71c0a6b98a87335b965015ef77c35 100644 (file)
 
 void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
 {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       lev =0; // get rid of warnings...
-#else
        FSGR_FORIJK_BOUNDS(lev) {
                if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
                        if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
@@ -50,15 +46,10 @@ void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
                }
        } // } TEST DEBUG
        if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
-#endif //! LBM_NOADCOARSENING==1
 }
        
 void LbmFsgrSolver::coarseAdvance(int lev)
 {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       lev =0; // get rid of warnings...
-#else
        LbmFloat calcCurrentMass = 0.0;
        LbmFloat calcCurrentVolume = 0.0;
 
@@ -177,10 +168,9 @@ void LbmFsgrSolver::coarseAdvance(int lev)
   mLevel[lev].lmass   = calcCurrentMass   * mLevel[lev].lcellfactor;
   mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
 #if ELBEEM_PLUGIN!=1
-  errMsg("DFINI", " m l"<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor );
-  errMsg("DFINI", " v l"<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor );
+  debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "mass: lev="<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
+  debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "volume: lev="<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
 #endif // ELBEEM_PLUGIN
-#endif //! LBM_NOADCOARSENING==1
 }
 
 /*****************************************************************************/
@@ -191,10 +181,6 @@ void LbmFsgrSolver::coarseAdvance(int lev)
 // get dfs from level (lev+1) to (lev) coarse border nodes
 void LbmFsgrSolver::coarseRestrictFromFine(int lev)
 {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       lev =0; // get rid of warnings...
-#else
        if((lev<0) || ((lev+1)>mMaxRefine)) return;
 #if FSGR_STRICT_DEBUG==1
        // reset all unused cell values to invalid
@@ -245,15 +231,9 @@ void LbmFsgrSolver::coarseRestrictFromFine(int lev)
                } // & fluid
        }}}
        if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
-#endif //! LBM_NOADCOARSENING==1
 }
 
 bool LbmFsgrSolver::adaptGrid(int lev) {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       lev =0; // get rid of warnings...
-       return false;
-#else
        if((lev<0) || ((lev+1)>mMaxRefine)) return false;
        bool change = false;
        { // refinement, PASS 1-3
@@ -706,7 +686,6 @@ bool LbmFsgrSolver::adaptGrid(int lev) {
 
        if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); }
        return change;
-#endif //! LBM_NOADCOARSENING==1
 }
 
 /*****************************************************************************/
@@ -715,10 +694,6 @@ bool LbmFsgrSolver::adaptGrid(int lev) {
 
 void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet)
 {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       i=j=k=srcSet=dstSet=lev =0; // get rid of warnings...
-#else
        LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
        LbmFloat *tcel = RACPNT(lev  , i,j,k      ,dstSet);
 
@@ -862,15 +837,9 @@ void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, i
        RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
        RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
 #                              endif // OPT3D==0
-#endif //! LBM_NOADCOARSENING==1
 }
 
 void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
-#if LBM_NOADCOARSENING==1
-       if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
-       i=j=k=dstSet=lev =0; // get rid of warnings...
-       t=0.0; flagSet=0; markNbs=false;
-#else
        LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
        LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
 
@@ -952,7 +921,6 @@ void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int d
 
        IDF_WRITEBACK;
        return;
-#endif //! LBM_NOADCOARSENING==1
 }
 
 
@@ -1008,8 +976,9 @@ void LbmFsgrSolver::adaptTimestep() {
 
        LbmFloat dtdiff = fabs(newdt - this->mpParam->getTimestep());
        if(!this->mSilent) {
-               debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt<<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff<<
-                       " simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
+               debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt
+                       <<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff
+                       <<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
 
        // in range, and more than X% change?
        //if( newdt <  this->mpParam->getTimestep() ) // DEBUG
@@ -1025,7 +994,8 @@ void LbmFsgrSolver::adaptTimestep() {
                        rescale = true;
                        if(!this->mSilent) {
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10);
-                               debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<this->mpParam->getTimestep()<<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
+                               debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<this->mpParam->getTimestep()
+                                               <<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<<
                                                "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
                        }
index 46462ec52deb05bc424a704aac6818b3e6808ccd..578e1079cd73d33e45eb006aaba3e1cb4c4b65fa 100644 (file)
 #endif
 #endif
 
+#if LBM_INCLUDE_TESTSOLVERS==1
+#include "solver_test.h"
+#endif // LBM_INCLUDE_TESTSOLVERS==1
+
 /*****************************************************************************/
 /*! cell access classes */
 class UniformFsgrCellIdentifier : 
@@ -216,10 +220,10 @@ class LbmFsgrSolver :
                //! notify object that dump is in progress (e.g. for field dump) 
                virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
 
-#if LBM_USE_GUI==1
+#              if LBM_USE_GUI==1
                //! show simulation info (implement LbmSolverInterface pure virtual func)
                virtual void debugDisplay(int set);
-#endif
+#              endif
                
                // implement CellIterator<UniformFsgrCellIdentifier> interface
                typedef UniformFsgrCellIdentifier stdCellId;
@@ -254,6 +258,7 @@ class LbmFsgrSolver :
                LBM_INLINED void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
                LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
                LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
+               LBM_INLINED void forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
                //! interpolate velocity and density at a given position
                void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz);
 
@@ -279,11 +284,11 @@ class LbmFsgrSolver :
                virtual vector<ntlGeometryObject*> getDebugObjects();
        
                // gui/output debugging functions
-#if LBM_USE_GUI==1
+#              if LBM_USE_GUI==1
                virtual void debugDisplayNode(int dispset, CellIdentifierInterface* cell );
                virtual void lbmDebugDisplay(int dispset);
                virtual void lbmMarkedCellDisplay();
-#endif // LBM_USE_GUI==1
+#              endif // LBM_USE_GUI==1
                virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1);
 
                //! for raytracing, preprocess
@@ -298,8 +303,10 @@ class LbmFsgrSolver :
                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);
-
+               
+               // loop over grid, stream&collide update
                void mainLoop(int lev);
+               // change time step size
                void adaptTimestep();
                //! init mObjectSpeeds for current parametrization
                void recalculateObjectSpeeds();
@@ -321,6 +328,11 @@ class LbmFsgrSolver :
                LBM_INLINED int getForZMaxBnd(int lev);
                LBM_INLINED int getForZMax1(int lev);
 
+               // touch grid and flags once
+               void preinitGrids();
+               // one relaxation step for standing fluid
+               void standingFluidPreinit();
+
 
                // member vars
 
@@ -355,6 +367,20 @@ class LbmFsgrSolver :
                LbmFloat mGfxEndTime;
                //! smoother surface initialization?
                int mInitSurfaceSmoothing;
+               //! surface generation settings, default is all off (=0)
+               //  each flag switches side on off,  fssgNoObs is for obstacle sides
+               //  -1 equals all on
+               typedef enum {
+                        fssgNormal   =  0,
+                        fssgNoNorth  =  1,
+                        fssgNoSouth  =  2,
+                        fssgNoEast   =  4,
+                        fssgNoWest   =  8,
+                        fssgNoTop    = 16,
+                        fssgNoBottom = 32,
+                        fssgNoObs    = 64
+               } fsSurfaceGen;
+               int mFsSurfGenSetting;
 
                //! lock time step down switching
                int mTimestepReduceLock;
@@ -392,10 +418,6 @@ class LbmFsgrSolver :
                int mMaxNoCells, mMinNoCells;
                LONGINT mAvgNumUsedCells;
 
-               //! 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
@@ -406,6 +428,7 @@ class LbmFsgrSolver :
                //! permanent movobj vert storage
          vector<ntlVec3Gfx>  mMOIVertices;
        vector<ntlVec3Gfx>  mMOIVerticesOld;
+         vector<ntlVec3Gfx>  mMOINormals;
 
                //! get isofield weights
                int mIsoWeightMethod;
@@ -443,6 +466,8 @@ class LbmFsgrSolver :
 
                //! debug function to disable standing f init
                int mDisableStandingFluidInit;
+               //! init 2d with skipped Y/Z coords
+               bool mInit2dYZ;
                //! debug function to force tadap syncing
                int mForceTadapRefine;
                //! border cutoff value
@@ -463,14 +488,31 @@ 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();
+               void cpDebugDisplay(int dispset);
+
+               void testXchng(); 
+       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
+               // former LbmModelLBGK  functions
                // relaxation funtions - implemented together with relax macros
                static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz);
                static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz);
                inline LbmFloat getLesNoneqTensorCoeff( LbmFloat df[],                          LbmFloat feq[] );
                inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo);
-               inline void collideArrays( int i, int j, int k, // position - more for debugging
+               inline void collideArrays( int lev, int i, int j, int k, // position - more for debugging
                                LbmFloat df[], LbmFloat &outrho, // out only!
                                // velocity modifiers (returns actual velocity!)
                                LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
@@ -479,12 +521,10 @@ class LbmFsgrSolver :
 
 
                // former LBM models
-       public:
-
-//! shorten static const definitions
-#define STCON static const
+               //! shorten static const definitions
+#              define STCON static const
 
-#if LBMDIM==3
+#              if LBMDIM==3
                
                //! id string of solver
                virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D3Q19]"); }
@@ -581,7 +621,7 @@ class LbmFsgrSolver :
                static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
                static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
 
-#else // end LBMDIM==3 , LBMDIM==2
+#              else // end LBMDIM==3 , LBMDIM==2
                
                //! id string of solver
                virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D2Q9]"); }
@@ -667,7 +707,7 @@ class LbmFsgrSolver :
                static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
                static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
 
-#endif  // LBMDIM==2
+#              endif  // LBMDIM==2
 };
 
 #undef STCON
@@ -691,7 +731,8 @@ class LbmFsgrSolver :
 
 // flag array defines -----------------------------------------------------------------------------------------------
 
-// lbm testsolver get index define
+// lbm testsolver get index define, note - ignores is (set) as flag
+// array is only a single entry
 #define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
 
 //! flag array acces macro
@@ -701,7 +742,6 @@ class LbmFsgrSolver :
 
 // array handling  -----------------------------------------------------------------------------------------------
 
-//#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)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set, l)*dTotalNum +(l)])
@@ -835,12 +875,17 @@ inline LbmFloat LbmFsgrSolver::getCollideEq(int l, LbmFloat rho,  LbmFloat ux, L
                                + rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
                                + 3.0 *tmp 
                                + 9.0/2.0 *(tmp*tmp) )
-                       );
+                       ); // */
 };
 
 /*****************************************************************************/
 /* init a given cell with flag, density, mass and equilibrium dist. funcs */
 
+void LbmFsgrSolver::forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
+       // also overwrite persisting flags
+       // function is useful for tracking accesses...
+       RFLAG(level,xx,yy,zz,set) = newflag;
+}
 void LbmFsgrSolver::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;
@@ -857,7 +902,6 @@ LbmFsgrSolver::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, Lb
        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;
@@ -875,7 +919,6 @@ LbmFsgrSolver::initVelocityCell(int level, int i,int j,int k, CellFlagType flag,
        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;
index 5fb3f84743a4fc4a479e144eea0a277a709e5922..1a19e1deabf0f9f8d03257ad9becfd2e55e084a2 100644 (file)
 // all boundary types at once
 #define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
 
+// helper for 2d init
+#define SWAPYZ(vec) { \
+               const LbmFloat tmp = (vec)[2]; \
+               (vec)[2] = (vec)[1]; (vec)[1] = tmp; }
+
 
 /*****************************************************************************/
 //! common variables 
@@ -302,15 +307,14 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mpPreviewSurface(NULL), 
        mTimeAdap(true), mForceTimeStepReduce(false),
        mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
-       mInitSurfaceSmoothing(0),
+       mInitSurfaceSmoothing(0), mFsSurfGenSetting(0),
        mTimestepReduceLock(0),
        mTimeSwitchCounts(0), mTimeMaxvelStepCnt(0),
        mSimulationTime(0.0), mLastSimTime(0.0),
        mMinTimestep(0.0), mMaxTimestep(0.0),
        mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
-       mDropMode(1), mDropSize(0.15), mDropSpeed(0.0),
        mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
-       mMOIVertices(), mMOIVerticesOld(),
+       mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
        mIsoWeightMethod(1),
        mMaxRefine(1), 
        mDfScaleUp(-1.0), mDfScaleDown(-1.0),
@@ -319,6 +323,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
        mLastOmega(1e10), mLastGravity(1e10),
        mNumInvIfTotal(0), mNumFsgrChanges(0),
        mDisableStandingFluidInit(0),
+       mInit2dYZ(false),
        mForceTadapRefine(-1), mCutoff(-1)
 {
   // not much to do here... 
@@ -456,12 +461,21 @@ void LbmFsgrSolver::parseAttrList()
        this->mSmoothSurface = this->mpAttrs->readFloat("smoothsurface", this->mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
        this->mSmoothNormals = this->mpAttrs->readFloat("smoothnormals", this->mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
 
+       mFsSurfGenSetting = this->mpAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
+       if(mFsSurfGenSetting==-1) {
+               // all on
+               mFsSurfGenSetting = 
+                        fssgNormal   | fssgNoNorth  | fssgNoSouth  | fssgNoEast   |
+                        fssgNoWest   | fssgNoTop    | fssgNoBottom | fssgNoObs   ;
+       }
+
        // refinement
        mMaxRefine = this->mRefinementDesired;
        mMaxRefine  = this->mpAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
        if(mMaxRefine<0) mMaxRefine=0;
        if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
        mDisableStandingFluidInit = this->mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
+       mInit2dYZ = this->mpAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
        mForceTadapRefine = this->mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
 
        // demo mode settings
@@ -489,7 +503,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 );
@@ -513,6 +527,7 @@ void LbmFsgrSolver::initLevelOmegas()
        this->mOmega = this->mpParam->calculateOmega(mSimulationTime);
        this->mGravity = vec2L( this->mpParam->calculateGravity(mSimulationTime) );
        this->mSurfaceTension = 0.; //this->mpParam->calculateSurfaceTension(); // unused
+       if(mInit2dYZ) { SWAPYZ(mGravity); }
 
        // check if last init was ok
        LbmFloat gravDelta = norm(this->mGravity-mLastGravity);
@@ -831,16 +846,33 @@ bool LbmFsgrSolver::initializeSolverMemory()
                isoend[2] = se;
                twodOff = 2;
        }
+       int isosx = this->mSizex+2;
+       int isosy = this->mSizey+2;
+       int isosz = this->mSizez+2+twodOff;
+
+       // MPT
+#if LBM_INCLUDE_TESTSOLVERS==1
+       if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
+               int xfac = 2;
+               isosx *= xfac;
+               isoend[0] = isostart[0] + (isoend[0]-isostart[0])*(LbmFloat)xfac;
+       }
+#endif // LBM_INCLUDE_TESTSOLVERS==1
+
        //errMsg(" SETISO ", " "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5)<<" "<<(LbmFloat)(this->mSizex+1.0)<<" " );
-       this->mpIso->setStart( vec2G(isostart) );
-       this->mpIso->setEnd(   vec2G(isoend) );
+       mpIso->setStart( vec2G(isostart) );
+       mpIso->setEnd(   vec2G(isoend) );
        LbmVec isodist = isoend-isostart;
-       this->mpIso->initializeIsosurface( this->mSizex+2, this->mSizey+2, this->mSizez+2+twodOff, vec2G(isodist) );
-       for(int ak=0;ak<this->mSizez+2+twodOff;ak++) 
-               for(int aj=0;aj<this->mSizey+2;aj++) 
-                       for(int ai=0;ai<this->mSizex+2;ai++) { *this->mpIso->getData(ai,aj,ak) = 0.0; }
+       mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
+
+       // reset iso field
+       for(int ak=0;ak<isosz;ak++) 
+               for(int aj=0;aj<isosy;aj++) 
+                       for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
+
 
   /* init array (set all invalid first) */
+       preinitGrids();
        for(int lev=0; lev<=mMaxRefine; lev++) {
                FSGR_FORIJK_BOUNDS(lev) {
                        RFLAG(lev,i,j,k,0) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
@@ -968,6 +1000,45 @@ bool LbmFsgrSolver::initializeSolverGrids() {
                        }
        }
        // Symmetry tests */
+       // vortt
+#if LBM_INCLUDE_TESTSOLVERS==1
+       if(( strstr( this->getName().c_str(), "vorttfluid" ) != NULL) && (LBMDIM==2)) {
+               errMsg("VORTT","init");
+               int level=mMaxRefine;
+               int cx = mLevel[level].lSizex/2;
+               int cyo = mLevel[level].lSizey/2;
+               int sx = mLevel[level].lSizex/8;
+               int sy = mLevel[level].lSizey/8;
+               LbmFloat rho = 1.;
+               LbmFloat rhomass = 1.;
+               LbmFloat uFactor = 0.15;
+               LbmFloat vdist = 1.0;
+
+               int cy1=cyo-(int)(vdist*sy);
+               int cy2=cyo+(int)(vdist*sy);
+
+               //for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) {      
+               for(int j=1;j<mLevel[level].lSizey-1;j++)
+                       for(int i=1;i<mLevel[level].lSizex-1;i++) {
+                               LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
+                               LbmFloat d2 = norm(LbmVec(cx,cy2,0.)-LbmVec(i,j,0));
+                               bool in1 = (d1<=(LbmFloat)(sx));
+                               bool in2 = (d2<=(LbmFloat)(sx));
+                               LbmVec uvec(0.);
+                         LbmVec v1 = getNormalized( cross( LbmVec(cx,cy1,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
+                         LbmVec v2 = getNormalized( cross( LbmVec(cx,cy2,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
+                               LbmFloat w1=1., w2=1.;
+                               if(!in1) w1=(LbmFloat)(sx)/(1.5*d1);
+                               if(!in2) w2=(LbmFloat)(sx)/(1.5*d2);
+                               if(!in1) w1=0.; if(!in2) w2=0.; // sharp falloff
+                         uvec += v1*w1;
+                         uvec += v2*w2;
+                               initVelocityCell(level, i,j,0, CFFluid, rho, rhomass, uvec );
+                               //errMsg("VORTT","init uvec"<<uvec);
+                       }
+
+       }
+#endif // LBM_INCLUDE_TESTSOLVERS==1
 
        //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
 
@@ -1083,13 +1154,14 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
 
 #define POS2GRID_CHECK(vec,n) \
                                monTotal++;\
+                               int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
+                               if(k!=0) continue; \
                                const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
                                if(i<=0) continue; \
                                if(i>=mLevel[level].lSizex-1) continue; \
                                const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
                                if(j<=0) continue; \
                                if(j>=mLevel[level].lSizey-1) continue;  \
-                               const int k=0; \
 
 #else // LBMDIM -> 3
 #define POS2GRID_CHECK(vec,n) \
@@ -1117,7 +1189,25 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
                                                if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
                                                if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
                                        } \
-                               } }
+                               } } \
+                               if(ntype&(CFBndFreeslip)) { \
+                                       const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
+                                       const LbmVec oldov=objvel; /*DEBUG*/ \
+                                       objvel = vec2L((*pNormals)[n]) *dp; \
+                                       /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
+                               } \
+                               else if(ntype&(CFBndPartslip)) { \
+                                       const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
+                                       const LbmVec oldov=objvel; /*DEBUG*/ \
+                                       /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
+                                       const LbmFloat partv = mObjectPartslips[OId]; \
+                                       /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, this->dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
+                                       /* m[l] = (RAC(ccel, this->dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
+                                       objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
+                               }
+
+#define TTT \
+
 
 /*****************************************************************************/
 //! init moving obstacles for next sim step sim 
@@ -1128,6 +1218,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        // movobj init
        const int level = mMaxRefine;
        const int workSet = mLevel[level].setCurr;
+       const int otherSet = mLevel[level].setOther;
        LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
        // for debugging - check targetTime check during DEFAULT STREAM
        LbmFloat targetTime = mSimulationTime + this->mpParam->getTimestep();
@@ -1180,6 +1271,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
 
                        otype = ntype = CFInvalid;
                        switch(obj->getGeoInitType()) {
+                                       /*
                                case FGI_BNDPART: 
                                case FGI_BNDFREE: 
                                        if(!staticInit) {
@@ -1191,16 +1283,15 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                        }
                                        break; 
                                        // off */
-                                       /*
                                case FGI_BNDPART: rhomass = BND_FILL;
-                                       otype = ntype = CFBnd|CFBndPartslip;
+                                       otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
                                        break;
                                case FGI_BNDFREE: rhomass = BND_FILL;
-                                       otype = ntype = CFBnd|CFBndFreeslip;
+                                       otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
                                        break;
                                        // off */
                                case FGI_BNDNO:   rhomass = BND_FILL;
-                                       otype = ntype = CFBnd|CFBndNoslip;
+                                       otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
                                        break;
                                case FGI_FLUID: 
                                        otype = ntype = CFFluid; 
@@ -1223,26 +1314,32 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                        mObjectSpeeds[OId] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
                        debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
 
+                       //vector<ntlVec3Gfx> tNormals;
+                       vector<ntlVec3Gfx> *pNormals = NULL;
+                       mMOINormals.clear();
+                       if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
+
                        mMOIVertices.clear();
                        if(obj->getMeshAnimated()) { 
                                // do two full update
+                               // TODO tNormals handling!?
                                mMOIVerticesOld.clear();
-                               obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices,  mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
+                               obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
                                monTrafo += mMOIVerticesOld.size();
-                               obj->applyTransformation(sourceTime, &mMOIVerticesOld,NULL, 0, mMOIVerticesOld.size(), false );
+                               obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
                                monTrafo += mMOIVertices.size();
-                               obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
+                               obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
                        } else {
                                // only do transform update
-                               obj->getMovingPoints(mMOIVertices);
+                               obj->getMovingPoints(mMOIVertices,pNormals);
                                mMOIVerticesOld = mMOIVertices;
                                // WARNING - assumes mSimulationTime is global!?
-                               obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
+                               obj->applyTransformation(targetTime, &mMOIVertices,pNormals, 0, mMOIVertices.size(), false );
                                monTrafo += mMOIVertices.size();
 
                                // correct flags from last position, but extrapolate
                                // velocity to next timestep
-                               obj->applyTransformation(sourceTime, &mMOIVerticesOld,NULL, 0, mMOIVerticesOld.size(), false );
+                               obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
                                monTrafo += mMOIVerticesOld.size();
                        }
 
@@ -1263,33 +1360,22 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                LbmFloat massCheck = 0.;
                                int massReinits=0;                              
                                bool fillCells = (mObjectMassMovnd[OId]<=-1.);
-                               LbmFloat massCScale; //, massCScalePos, massCScaleNeg;
-                               //if(mInitialMass>0.) {
-                                       //massCScale = (mInitialMass-mObjectMassMovnd[OId])/mInitialMass;
-                                       massCScale = 1.-(mObjectMassMovnd[OId]/(LbmFloat)(mMOIVertices.size()/10) );
-                                       //massCScalePos = MIN(massCScale,1.);
-                                       //massCScaleNeg = MAX(massCScale,1.);
-                               //} else {
-                                       //massCScale = massCScalePos = massCScaleNeg = 1.;
-                               //}
 
                                // first pass - set new obs. cells
                                if(active) {
                                        for(size_t n=0; n<mMOIVertices.size(); n++) {
-                                               //errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
+                                               //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
                                                POS2GRID_CHECK(mMOIVertices,n);
-                                               //if(i==30 && j==14) { errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<" "); }
+                                               //{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
                                                if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
                                                monPoints++;
                                                
                                                // check mass
                                                if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
-                                                       FORDF0 {
-                                                               massCheck -= QCELL(level, i,j,k, workSet, l);
-                                                       }
+                                                       FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
                                                        massReinits++;
                                                }
-                                               if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
+                                               else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
                                                        massCheck -= QCELL(level, i,j,k, workSet, dMass);
                                                        massReinits++;
                                                }
@@ -1314,60 +1400,94 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                                        const LbmFloat ux = this->dfDvecX[l]*objvel[0];
                                                                        const LbmFloat uy = this->dfDvecY[l]*objvel[1];
                                                                        const LbmFloat uz = this->dfDvecZ[l]*objvel[2];
-                                                                       LbmFloat factor = 2.0*this->dfLength[l]* 3.0 *(ux+uy+uz); // rhoTest, dont multiply by density...
+
+                                                                       /*LbmFloat *rhodstCell = RACPNT(level, i+dfVecX[l],j+dfVecY[l],k+dfVecZ[l],workSet);
+                                                                       LbmFloat rhodst = RAC(rhodstCell,dC);
+                                                                       for(int ll=1; ll<19; ll++) { rhodst += RAC(rhodstCell,ll); }
+                                                                       rhodst = 1.; // DEBUG TEST! */
+
+                                                                       LbmFloat factor = 2. * this->dfLength[l] * 3.0 * (ux+uy+uz); // 
+                                                                               /* FSTEST 
+                                                                       */
+                                                                       if(ntype&(CFBndFreeslip|CFBndPartslip)) {
+                                                                               // missing, diag mass checks...
+                                                                               //if(l<=LBMDIM*2) massCheck += factor;
+
+                                                                               //if(l<=LBMDIM*2) factor *= 1.4142;
+                                                                               //if(l>LBMDIM*2) factor = 0.;
+                                                                               //else 
+                                                                               //factor = 0.; // TODO, optimize
+                                                                               factor *= 2.0; // TODO, optimize
+                                                                       } else {
+                                                                               factor *= 1.2; // TODO, optimize
+                                                                       }
                                                                        RAC(dstCell,l) = factor;
-                                                                       massCheck += RAC(dstCell,l); 
+                                                                       massCheck += factor;
                                                                } else {
                                                                        RAC(dstCell,l) = 0.;
                                                                }
                                                        }
+
+#if NEWDIRVELMOTEST==1
+                                                       FORDF1 { RAC(dstCell,l) = 0.; }
+                                                       RAC(dstCell,dMass)  = objvel[0];
+                                                       RAC(dstCell,dFfrac) = objvel[1];
+                                                       RAC(dstCell,dC)     = objvel[2];
+#endif // NEWDIRVELMOTEST==1
                                                } else {
                                                        FORDF1 { RAC(dstCell,l) = 0.0; }
                                                }
-                                               //errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" objvel"<<objvel<<" ul"<<PRINT_VEC(ux,uy,uz) );
                                                RAC(dstCell, dFlux) = targetTime;
+                                               //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" dflt="<<RAC(dstCell, dFlux) );
                                                monObsts++;
                                        } // points
                                } // bnd, is active?
 
                                // second pass, remove old ones
+                               // warning - initEmptyCell et. al dont overwrite OId or persist flags...
                                if(wasActive) {
                                        for(size_t n=0; n<mMOIVerticesOld.size(); n++) {
                                                POS2GRID_CHECK(mMOIVerticesOld,n);
                                                monPoints++;
                                                if((RFLAG(level, i,j,k, workSet) == otype) &&
                                                         (QCELL(level, i,j,k, workSet, dFlux) != targetTime)) {
-                                                       //? unused ntlVec3Gfx objvel= (mMOIVertices[n]-mMOIVerticesOld[n]);
                                                        // from mainloop
                                                        nbored = 0;
+                                                       // TODO: optimize for OPT3D==0
                                                        FORDF1 {
-                                                               rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
+                                                               //rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
+                                                               rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance
                                                                nbored |= rflagnb[l];
                                                        } 
                                                        CellFlagType settype = CFInvalid;
-                                                       //LbmFloat avgrho=0.0, avgux=0.0, avguy=0.0, avguz=0.0; 
                                                        if(nbored&CFFluid) {
+                                                       //if(nbored&(CFFluid|CFInter)) {
                                                                settype = CFInter|CFNoInterpolSrc; 
-                                                               if(fillCells) rhomass = 1.0;
-                                                               else rhomass = 0.;
+                                                               rhomass = 1.5;
+                                                               if(!fillCells) rhomass = 0.;
+                                                               //settype = CFFluid|CFNoInterpolSrc; rhomass=1.; // test
+                                                               //rhomass = 1.01; // DEBUGT
 
-                                                               //interpolateCellValues(level,i,j,k, workSet, avgrho,avgux,avguy,avguz);
-                                                               //LbmVec speed(avgux,avguy,avguz);
-                                                               //initVelocityCell(level, i,j,k, settype, avgrho, rhomass, speed );
                                                                OBJVEL_CALC;
-                                                               initVelocityCell(level, i,j,k, settype, 1., rhomass, objvel );
+                                                               if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
+                                                               //settype=CFFluid; // rhomass = 1.; objvel = LbmVec(0.); // DEBUG TEST
+
+                                                               // new interpolate values
+                                                               LbmFloat avgrho = 0.0;
+                                                               LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
+                                                               interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
+                                                               //objvel = LbmVec(avgux,avguy,avguz); 
+                                                               //avgrho=1.;
+                                                               initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
+                                                               //errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
+
+                                                               // old - fixed init
+                                                               //initVelocityCell(level, i,j,k, settype, 1., rhomass, objvel );
                                                                massCheck += rhomass;
-                                                       } 
-                                                       /*else if((nbored&CFInter)&&(fillCells)) {
-                                                               settype = CFInter|CFNoInterpolSrc; rhomass = 1.0;
-                                                               _interpolateCellValues(level,i,j,k, workSet, avgrho,avgux,avguy,avguz);
-                                                       }  // */
-                                                       else {
+                                                       } else {
                                                                settype = CFEmpty; rhomass = 0.0;
                                                                initEmptyCell(level, i,j,k, settype, 1., rhomass );
                                                        }
-                                                       //settype = CFBnd|CFBndNoslip; rhomass = 0.0;
-                                                       //avgux=avguy=avguz=0.0; avgrho=1.0;
                                                        monFluids++;
                                                        massReinits++;
                                                } // flag & simtime
@@ -1376,8 +1496,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
 
                                // only compute mass transfer when init is done
                                if(this->mInitDone) {
-                                       errMsg("initMov","Massd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
-                                               " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" massCScale"<<massCScale<<" inim:"<<mInitialMass ); 
+                                       errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
+                                               " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass ); 
                                        mObjectMassMovnd[OId] += massCheck;
                                }
                        } // bnd, active
@@ -1437,9 +1557,10 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                POS2GRID_CHECK(mMOIVertices,n);
                                                if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
                                                if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
-                                                       changeFlag(level, i,j,k, workSet, set2Flag);
+                                                       forceChangeFlag(level, i,j,k, workSet, set2Flag);
                                                } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
-                                                       changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag);
+                                                       forceChangeFlag(level, i,j,k, workSet, 
+                                                                       (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
                                                }
                                        }
                                } // second static inflow pass
@@ -1453,7 +1574,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                        // FIXME check fluid/inter cells for non-static!?
                                        if(!(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter))) {
                                                if((staticInit)&&(RFLAG(level, i,j,k, workSet)==CFEmpty)) {
-                                                       changeFlag(level, i,j,k, workSet, CFMbndOutflow); }
+                                                       forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
                                                continue;
                                        }
                                        monFluids++;
@@ -1481,9 +1602,10 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
                                                POS2GRID_CHECK(mMOIVertices,n);
                                                if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
                                                if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
-                                                       changeFlag(level, i,j,k, workSet, set2Flag);
+                                                       forceChangeFlag(level, i,j,k, workSet, set2Flag);
                                                } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
-                                                       changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag);
+                                                       forceChangeFlag(level, i,j,k, workSet, 
+                                                                       (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
                                                }
                                        }
                                } // second static outflow pass
@@ -1511,6 +1633,15 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
 }
 
 
+// geoinit position
+
+#define GETPOS(i,j,k)  \
+       ntlVec3Gfx ggpos = \
+               ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
+                                 iniPos[1]+ dvec[1]*(gfxReal)(j), \
+                                 iniPos[2]+ dvec[2]*(gfxReal)(k) ); \
+  if((LBMDIM==2)&&(mInit2dYZ)) { SWAPYZ(ggpos); }
+
 /*****************************************************************************/
 /*! perform geometry init (if switched on) */
 /*****************************************************************************/
@@ -1580,6 +1711,7 @@ bool LbmFsgrSolver::initGeometryFlags() {
        if(LBMDIM==2) {
                dvec[2] = 0.0; 
                iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))+(dvec*0.5);
+               //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
        } else {
                iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
        }
@@ -1587,16 +1719,13 @@ bool LbmFsgrSolver::initGeometryFlags() {
 
        // first init boundary conditions
        // invalid cells are set to empty afterwards
-#define GETPOS(i,j,k) \
-                                               ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
-                                                           iniPos[1]+ dvec[1]*(gfxReal)(j), \
-                                                           iniPos[2]+ dvec[2]*(gfxReal)(k) )
        for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
                for(int j=1;j<mLevel[level].lSizey-1;j++) {
                        for(int i=1;i<mLevel[level].lSizex-1;i++) {
                                ntype = CFInvalid;
                                
-                               const bool inside = this->geoInitCheckPointInside( GETPOS(i,j,k) , FGI_ALLBOUNDS, OId, distance);
+                               GETPOS(i,j,k);
+                               const bool inside = this->geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
                                if(inside) {
                                        pObj = (*this->mpGiObjects)[OId];
                                        switch( pObj->getGeoInitType() ){
@@ -1620,6 +1749,9 @@ bool LbmFsgrSolver::initGeometryFlags() {
                                                rhomass = BND_FILL;
                                                ntype = CFBnd|CFBndNoslip; 
                                                break;
+                                       case FGI_BNDPART: 
+                                               rhomass = BND_FILL;
+                                               ntype = CFBnd|CFBndPartslip; break;
                                        case FGI_BNDFREE: 
                                                rhomass = BND_FILL;
                                                ntype = CFBnd|CFBndFreeslip; break;
@@ -1662,9 +1794,8 @@ bool LbmFsgrSolver::initGeometryFlags() {
                                if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
                                ntype = CFInvalid;
                                int inits = 0;
-                               //if((i==1) && (j==31) && (k==48)) globGeoInitDebug=1;
-                               //else globGeoInitDebug=0;
-                               const bool inside = this->geoInitCheckPointInside( GETPOS(i,j,k) , FGI_FLUID, OId, distance);
+                               GETPOS(i,j,k);
+                               const bool inside = this->geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
                                if(inside) {
                                        ntype = CFFluid;
                                }
@@ -1723,23 +1854,16 @@ void LbmFsgrSolver::initFreeSurfaces() {
        // set interface cells 
        FSGR_FORIJK1(mMaxRefine) {
                if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid )) {
-                       //int initInter = 0; // check for neighboring empty cells 
                        FORDF1 {
                                int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
                                if( TESTFLAG( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr), CFEmpty ) ) {
                                        LbmFloat arho=0., aux=0., auy=0., auz=0.;
                                        interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
                                        //errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
-                                       initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
+                                       // unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
                                        initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
-                                       //initEmptyCell(level, i,j,k, settype, avgrho, rhomass, speed ); */
-                                       //initInter = 1;
                                }
                        }
-                       /*if(initInter) {
-                               QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr, dMass) = interfaceFill;
-                               RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) = RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther) = CFInter;
-                       } // */
                }
        }
 
@@ -1950,10 +2074,6 @@ void LbmFsgrSolver::initStandingFluidGradient() {
        if(haveStandingFluid) {
                int rhoworkSet = mLevel[lev].setCurr;
                myTime_t timestart = getTime(); // FIXME use user time here?
-               LbmFloat lcsmqo;
-#if OPT3D==1 
-               LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
-#endif // OPT3D==true 
 
                GRAVLOOP {
                        int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
@@ -1980,87 +2100,21 @@ void LbmFsgrSolver::initStandingFluidGradient() {
                        }
 
                } // GRAVLOOP
+
                debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<<
                        (1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8);
                
-#if ELBEEM_PLUGIN!=1 && LBMDIM==3
-               /*int lowj = 0;
-               for(int k=1;k<mLevel[lev].lSizez-1;++k) {
-               for(int i=1;i<mLevel[lev].lSizex-1;++i) {
-                       LbmFloat rho = 1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
-                       RFLAG(lev, i,lowj,k, rhoworkSet^1) =
-                       RFLAG(lev, i,lowj,k, rhoworkSet) = CFFluid;
-                       FORDF0 { QCELL(lev, i,lowj,k, rhoworkSet, l) = this->dfEquil[l]*rho; }
-                       QCELL(lev, i,lowj,k, rhoworkSet, dMass) = rho;
-               } } // */
-#endif 
-
                int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
                preinitSteps = (haveStandingFluid>>2); // not much use...?
-               //preinitSteps = 4; // DEBUG!!!!
-               //this->mInitDone = 1; // GRAVTEST
                //preinitSteps = 0;
-               debMsgNnl("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
+               debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
                for(int s=0; s<preinitSteps; s++) {
-                       int workSet = SRCS(lev); //mLevel[lev].setCurr;
-                       int otherSet = TSET(lev); //mLevel[lev].setOther;
-                       debMsgDirect(".");
-                       if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "s="<<s<<" curset="<<workSet<<" srcs"<<SRCS(lev), 10);
-                       LbmFloat *ccel;
-                       LbmFloat *tcel;
-                       LbmFloat m[LBM_DFNUM];
-
-               // grav loop not necessary here
-#define NBFLAG(l) (nbflag[(l)])
-               LbmFloat rho, ux,uy,uz, usqr; 
-               int kstart=getForZMinBnd(), kend=getForZMaxBnd(mMaxRefine);
-#if COMPRESSGRIDS==0
-               for(int k=kstart;k<kend;++k) {
-#else // COMPRESSGRIDS==0
-               int kdir = 1; // COMPRT ON
-               if(mLevel[lev].setCurr==1) {
-                       kdir = -1;
-                       int temp = kend;
-                       kend = kstart-1;
-                       kstart = temp-1;
-               } // COMPRT
-               for(int k=kstart;k!=kend;k+=kdir) {
-#endif // COMPRESSGRIDS==0
-
-               for(int j=0;j<mLevel[lev].lSizey-0;++j) {
-               for(int i=0;i<mLevel[lev].lSizex-0;++i) {
-                               const CellFlagType currFlag = RFLAG(lev, i,j,k,workSet);
-                               if( (currFlag & (CFEmpty|CFBnd)) ) continue;
-                               ccel = RACPNT(lev, i,j,k,workSet); 
-                               tcel = RACPNT(lev, i,j,k,otherSet);
-
-                               if( (currFlag & (CFInter)) ) {
-                                       // copy all values
-                                       for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
-                                       continue;
-                               }
-
-                               if( (currFlag & CFNoBndFluid)) {
-                                       OPTIMIZED_STREAMCOLLIDE;
-                               } else {
-                                       FORDF1 {
-                                               nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
-                                       } 
-                                       DEFAULT_STREAM;
-                                       //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
-                                       DEFAULT_COLLIDEG(mLevel[lev].gravity);
-                               }
-                               for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
-                       } } } // GRAVLOOP
-
-                       mLevel[lev].setOther = mLevel[lev].setCurr;
-                       mLevel[lev].setCurr ^= 1;
+                       // in solver main cpp
+                       standingFluidPreinit();
                }
-               //this->mInitDone = 0;  // GRAVTEST
-               // */
 
                myTime_t timeend = getTime();
-               debMsgDirect(" done, "<<getTimeString(timeend-timestart)<<" \n");
+               debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
 #undef  NBFLAG
        }
 }
@@ -2150,17 +2204,15 @@ LbmFsgrSolver::interpolateCellValues(
        LbmFloat avgrho = 0.0;
        LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
        LbmFloat cellcnt = 0.0;
-       //LbmFloat avgnbdf[LBM_DFNUM];
-       //FORDF0M { avgnbdf[m]= 0.0; }
 
        for(int nbl=1; nbl< this->cDfNum ; ++nbl) {
-               if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFFluid) || 
-                               ((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
-                                (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) )) { 
+               if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
+               if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){
+                               //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
+                                //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) { 
                        cellcnt += 1.0;
                        for(int rl=0; rl< this->cDfNum ; ++rl) { 
                                LbmFloat nbdf =  QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
-                               //avgnbdf[rl] += nbdf;
                                avgux  += (this->dfDvecX[rl]*nbdf); 
                                avguy  += (this->dfDvecY[rl]*nbdf);  
                                avguz  += (this->dfDvecZ[rl]*nbdf);  
@@ -2171,20 +2223,18 @@ LbmFsgrSolver::interpolateCellValues(
 
        if(cellcnt<=0.0) {
                // no nbs? just use eq.
-               //FORDF0 { QCELL(level,ei,ej,ek, workSet, l) = this->dfEquil[l]; }
                avgrho = 1.0;
                avgux = avguy = avguz = 0.0;
                //TTT mNumProblems++;
 #if ELBEEM_PLUGIN!=1
                //this->mPanic=1; 
                // this can happen for worst case moving obj scenarios...
-               errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek)); //,SIMWORLD_GENERICERROR);
+               errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
 #endif // ELBEEM_PLUGIN
        } else {
                // init speed
                avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
                avgrho /= cellcnt;
-               //FORDF0M { avgnbdf[m] /= cellcnt; } // CHECK FIXME test?
        }
 
        retrho = avgrho;
index 7b3f71e438ceaf087298f347fac2586486134cba..e9db2a10db750976ac3481b0b44aca09e0d6f198 100644 (file)
@@ -285,7 +285,11 @@ void LbmSolverInterface::initGeoTree() {
        if(mpGiTree != NULL) delete mpGiTree;
        char treeFlag = (1<<(this->mLbmInitId+4));
        mpGiTree = new ntlTree( 
+# if FSGR_STRICT_DEBUG!=1
                        15, 8,  // TREEwarning - fixed values for depth & maxtriangles here...
+# else // FSGR_STRICT_DEBUG==1
+                       10, 20,  // reduced/slower debugging values
+# endif // FSGR_STRICT_DEBUG==1
                        scene, treeFlag );
 }
 
index 4cbdd2945e58fb683e5753ea1f74e021c254e79d..b95ce0358cbea95a8686b48140ce814814bbcb1e 100644 (file)
@@ -128,6 +128,7 @@ typedef int BubbleId;
 
 // above 24 is used to encode in/outflow object type
 #define CFPersistMask (0xFF000000 | CFMbndInflow | CFMbndOutflow)
+#define CFNoPersistMask (~CFPersistMask)
 
 
 // nk
index 8fa3b8d8a65456ba4e7cd12f80c57c82a1ccee27..964192420e4bbd4c50f715e9f75080fdd3f244c0 100644 (file)
 #include "solver_class.h"
 #include "solver_relax.h"
 #include "particletracer.h"
+#include "loop_tools.h"
 
 /*****************************************************************************/
 /*! perform a single LBM step */
 /*****************************************************************************/
 
+double globdfcnt;
+double globdfavg[19];
+double globdfmax[19];
+double globdfmin[19];
 
 void LbmFsgrSolver::step() { 
        initLevelOmegas();
        stepMain(); 
+
+       // intlbm test
+       if(0) {
+               if(this->mStepCnt<5) {
+                       // init
+                       globdfcnt=0.;
+                       FORDF0{ 
+                               globdfavg[l] = 0.;
+                               globdfmax[l] = -1000.; //this->dfEquil[l];
+                               globdfmin[l] = 1000.; // this->dfEquil[l];
+                       }
+               } else {
+
+                       int workSet = mLevel[mMaxRefine].setCurr;
+                       for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
+                               for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
+                                       for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
+                                               //float val = 0.;
+                                               if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) {
+                                                       //val = QCELL(mMaxRefine,i,j,k, workSet,dFfrac);
+                                                       FORDF0{ 
+                                                               const double df = (double)QCELL(mMaxRefine,i,j,k, workSet,l);
+                                                               globdfavg[l] += df;
+                                                               if(df>globdfmax[l]) globdfmax[l] = df;
+                                                               //if(df<=0.01) { errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l); }
+                                                               if(df<globdfmin[l]) globdfmin[l] = df;
+                                                               //errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l<<" "<<df<<" "<<globdfmin[l]); 
+                                                       }
+                                                       globdfcnt+=1.;
+                                               }
+                                       }
+                               }
+                       }
+                       if(this->mStepCnt%10==0) {
+                               FORDF0{ errMsg("GLOBDF","l="<<l<<" avg="<<(globdfavg[l]/globdfcnt)<<"   max="<<globdfmax[l]<<"   min="<<globdfmin[l]<<"   "<<globdfcnt); }
+                       }
+
+               }
+       }
+       // intlbm test */
 }
 
 void LbmFsgrSolver::stepMain()
@@ -189,6 +234,7 @@ void LbmFsgrSolver::stepMain()
 
 #if LBM_INCLUDE_TESTSOLVERS==1
        if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); }
+       testXchng();
 #endif
 
        // advance positions with current grid
@@ -213,8 +259,13 @@ void LbmFsgrSolver::stepMain()
 
        // output total step time
        timeend = getTime();
+       double mdelta = (lastMass-mCurrentMass);
+       if(ABS(mdelta)<1e-12) mdelta=0.;
        debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt
-                       <<": dccd="<< mCurrentMass<<",d"<<(lastMass-mCurrentMass)<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
+                       <<": dccd="<< mCurrentMass
+                       <<",d"<<mdelta
+                       <<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
+                       <<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
                        <<" totst:"<<getTimeString(timeend-timestart), 3);
        // nicer output
        debMsgDirect(std::endl); 
@@ -248,7 +299,7 @@ void LbmFsgrSolver::fineAdvance()
        // time adaptivity
        this->mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
        //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
-       if(!this->mSilent){ errMsg("fineAdvance"," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps) ); }
+       if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
 
        // update other set
   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
@@ -257,62 +308,117 @@ void LbmFsgrSolver::fineAdvance()
 
        // flag init... (work on current set, to simplify flag checks)
        reinitFlags( mLevel[mMaxRefine].setCurr );
-       if(!this->mSilent){ errMsg("fineAdvance"," flags reinit on set "<< mLevel[mMaxRefine].setCurr ); }
+       if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); }
+
+       // DEBUG VEL CHECK
+       if(0) {
+               int lev = mMaxRefine;
+               int workSet = mLevel[lev].setCurr;
+               int mi=0,mj=0,mk=0;
+               LbmFloat compRho=0.;
+               LbmFloat compUx=0.;
+               LbmFloat compUy=0.;
+               LbmFloat compUz=0.;
+               LbmFloat maxUlen=0.;
+               LbmVec maxU(0.);
+               LbmFloat maxRho=-100.;
+               int ri=0,rj=0,rk=0;
+
+               FSGR_FORIJK1(lev) {
+                       if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) {
+                               compRho=QCELL(lev, i,j,k,workSet, dC);
+                               compUx = compUy = compUz = 0.0;
+                               for(int l=1; l<this->cDfNum; l++) {
+                                       LbmFloat df = QCELL(lev, i,j,k,workSet, l);
+                                       compRho += df;
+                                       compUx  += (this->dfDvecX[l]*df);
+                                       compUy  += (this->dfDvecY[l]*df); 
+                                       compUz  += (this->dfDvecZ[l]*df); 
+                               } 
+                               LbmVec u(compUx,compUy,compUz);
+                               LbmFloat nu = norm(u);
+                               if(nu>maxUlen) {
+                                       maxU = u;
+                                       maxUlen = nu;
+                                       mi=i; mj=j; mk=k;
+                               }
+                               if(compRho>maxRho) {
+                                       maxRho=compRho;
+                                       ri=i; rj=j; rk=k;
+                               }
+                       } else {
+                               continue;
+                       }
+               }
+
+               errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU);
+               errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho);
+               printLbmCell(lev, 30,36,23, -1);
+       } // DEBUG VEL CHECK
+
 }
 
 
-/*****************************************************************************/
-//! fine step function
-/*****************************************************************************/
 
+// fine step defines
 
 // access to own dfs during step (may be changed to local array)
 #define MYDF(l) RAC(ccel, l)
 
+// drop model definitions
+#define RWVEL_THRESH 1.5
+#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
+
+#if LBMDIM==3
+                       // normal
+#define SLOWDOWNREGION (this->mSizez/4)
+#else // LBMDIM==2
+                       // off
+#define SLOWDOWNREGION 10 
+#endif // LBMDIM==2
+#define P_LCSMQO 0.01
+
+/*****************************************************************************/
+//! fine step function
+/*****************************************************************************/
 void 
 LbmFsgrSolver::mainLoop(int lev)
 {
        // loops over _only inner_ cells  -----------------------------------------------------------------------------------
-       LbmFloat calcCurrentMass = 0.0;
-       LbmFloat calcCurrentVolume = 0.0;
-       int      calcCellsFilled = this->mNumFilledCells;
-       int      calcCellsEmptied = this->mNumEmptiedCells;
-       int      calcNumUsedCells = this->mNumUsedCells;
+       
+       // slow surf regions smooth (if below)
+       const LbmFloat smoothStrength = 0.0; //0.01;
+       const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
+       const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
+
        const int cutMin  = 1;
        const int cutConst = mCutoff+2;
 
+
 #      if LBM_INCLUDE_TESTSOLVERS==1
        // 3d region off... quit
        if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; }
-#endif // ELBEEM_PLUGIN!=1
-       //printLbmCell(lev, 6,6,16, mLevel[lev].setCurr ); // DEBUG
+#      endif // ELBEEM_PLUGIN!=1
        
+  // main loop region
+       const bool doReduce = true;
+       const int gridLoopBound=1;
+       GRID_REGION_INIT();
 #if PARALLEL==1
-#include "paraloop.h"
+#include "paraloopstart.h"
+       GRID_REGION_START();
 #else // PARALLEL==1
-  { // main loop region
-       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 );
+       GRID_REGION_START();
 #endif // PARALLEL==1
 
-       // local to loop
+       // local to main
        CellFlagType nbflag[LBM_DFNUM]; 
-       LbmFloat *ccel = NULL;
-       LbmFloat *tcel = NULL;
-       int oldFlag;
-       int newFlag;
-       int nbored;
+       int oldFlag, newFlag, nbored;
        LbmFloat m[LBM_DFNUM];
        LbmFloat rho, ux, uy, uz, tmp, usqr;
-       LbmFloat mass, change, lcsmqo=0.0;
-       usqr = tmp = 0.0; 
-#if OPT3D==1 
-       LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
-#endif // OPT3D==true 
 
+       // smago vars
+       LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
        
        // ifempty cell conversion flags
        bool iffilled, ifemptied;
@@ -320,72 +426,22 @@ LbmFsgrSolver::mainLoop(int lev)
        int recons[LBM_DFNUM];   // reconstruct this DF?
        int numRecons;           // how many are reconstructed?
 
-       // slow surf regions smooth (if below)
-       const LbmFloat smoothStrength = 0.0; //0.01;
-       const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
-       const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
-       
-       CellFlagType *pFlagSrc;
-       CellFlagType *pFlagDst;
-       pFlagSrc = &RFLAG(lev, 0,1, kstart,SRCS(lev)); // omp
-       pFlagDst = &RFLAG(lev, 0,1, kstart,TSET(lev)); // omp
-       ccel = RACPNT(lev, 0,1, kstart ,SRCS(lev)); // omp
-       tcel = RACPNT(lev, 0,1, kstart ,TSET(lev)); // omp
-       //CellFlagType *pFlagTar = NULL;
-       int pFlagTarOff;
-       if(mLevel[lev].setOther==1) pFlagTarOff = mLevel[lev].lOffsz;
-       else pFlagTarOff = -mLevel[lev].lOffsz;
-#define ADVANCE_POINTERS(p)    \
-       ccel += (QCELLSTEP*(p));        \
-       tcel += (QCELLSTEP*(p));        \
-       pFlagSrc+= (p); \
-       pFlagDst+= (p); \
-       i+= (p);
+       LbmFloat mass=0., change=0., lcsmqo=0.;
+       rho= ux= uy= uz= usqr= tmp= 0.; 
+       lcsmqadd = lcsmomega = 0.;
+       FORDF0{ lcsmeq[l] = 0.; }
 
        // ---
        // now stream etc.
-
-       //{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIM<<" mlsz:"<<mLevel[mMaxRefine].lSizez ); } // DEBUG
-
        // use //template functions for 2D/3D
-#if COMPRESSGRIDS==0
-  for(int k=kstart;k<kend;++k) {
-  for(int j=1;j<mLevel[lev].lSizey-1;++j) {
-  for(int i=0;i<mLevel[lev].lSizex-2;   ) {
-#else // COMPRESSGRIDS==0
-       int kdir = 1; // COMPRT ON
-       if(mLevel[mMaxRefine].setCurr==1) {
-               kdir = -1;
-               int temp = kend;
-               kend = kstart-1;
-               kstart = temp-1;
-       } // COMPRT
-
-#if PARALLEL==0
-  //const int id = 0, Nthrds = 1;
-       const int jstart = 0;
-       const int jend   = mLevel[mMaxRefine].lSizey;
-//#endif // PARALLEL==1
-#else // PARALLEL==1
-       PARA_INITIALIZE();
-       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
-  pFlagSrc = &RFLAG(lev, 0, jstart, k, SRCS(lev)); // omp test // COMPRT ON
-  pFlagDst = &RFLAG(lev, 0, jstart, k, TSET(lev)); // omp test
-  ccel = RACPNT(lev,     0, jstart, k, SRCS(lev)); // omp test
-  tcel = RACPNT(lev,     0, jstart, k, TSET(lev)); // omp test // COMPRT ON
-
-  //for(int j=1;j<mLevel[lev].lSizey-1;++j) {
-  for(int j=jstart;j!=jend;++j) {
-  for(int i=0;i<mLevel[lev].lSizex-2;   ) {
-#endif // COMPRESSGRIDS==0
+       GRID_LOOP_START();
+               // loop start
+               // stream from current set to other, then collide and store
+               //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
 
-               ADVANCE_POINTERS(1); 
-#if FSGR_STRICT_DEBUG==1
+#              if FSGR_STRICT_DEBUG==1
+               // safety pointer check
                rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
                if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
                    (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
@@ -405,13 +461,11 @@ LbmFsgrSolver::mainLoop(int lev)
                                        ); 
                        CAUSE_PANIC;
                }       
-#endif
+#              endif
                oldFlag = *pFlagSrc;
-               //DEBUG if(LBMDIM==2) { errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<< RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<" "); }
-               // stream from current set to other, then collide and store
                
                // old INTCFCOARSETEST==1
-               if( (oldFlag & (CFGrFromCoarse)) ) {  // interpolateFineFromCoarse test!
+               if( (oldFlag & (CFGrFromCoarse)) ) { 
                        if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
                                FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
                        } else {
@@ -419,7 +473,7 @@ LbmFsgrSolver::mainLoop(int lev)
                                calcNumUsedCells++;
                        }
                        continue; // interpolateFineFromCoarse test!
-               } // interpolateFineFromCoarse test!  old INTCFCOARSETEST==1
+               } // interpolateFineFromCoarse test! 
        
                if(oldFlag & (CFMbndInflow)) {
                        // fluid & if are ok, fill if later on
@@ -434,9 +488,10 @@ LbmFsgrSolver::mainLoop(int lev)
                                RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
                                RAC(tcel, dFlux) = FLUX_INIT;
                                changeFlag(lev, i,j,k, TSET(lev), CFInter);
-                               calcCurrentMass += iniRho; calcCurrentVolume += 1.0; calcNumUsedCells++;
+                               calcCurrentMass += iniRho; 
+                               calcCurrentVolume += 1.0; 
+                               calcNumUsedCells++;
                                mInitialMass += iniRho;
-                               //errMsg("INFLOW_DEBUG","ini at "<<PRINT_IJK<<" v="<<vel<<" inirho="<<iniRho);
                                // dont treat cell until next step
                                continue;
                        } 
@@ -456,11 +511,6 @@ LbmFsgrSolver::mainLoop(int lev)
                                // same as ifemptied for if below
                                LbmPoint oemptyp; oemptyp.flag = 0;
                                oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
-#if PARALLEL==1
-                               //calcListEmpty[id].push_back( oemptyp );
-#else // PARALLEL==1
-                               //mListEmpty.push_back( oemptyp );
-#endif // PARALLEL==1
                                LIST_EMPTY(oemptyp);
                                calcCellsEmptied++;
                                continue;
@@ -597,7 +647,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
+               //onlyBndnb = false; // DEBUG test off
 
                FORDF1 { // dfl loop
                        recons[l] = 0;
@@ -690,13 +740,13 @@ LbmFsgrSolver::mainLoop(int lev)
                if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
                if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
                ny = 0.5* (nv2-nv1);
-#if LBMDIM==3
+#              if LBMDIM==3
                if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
                if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
                nz = 0.5* (nv2-nv1);
-#else // LBMDIM==3
+#              else // LBMDIM==3
                nz = 0.0;
-#endif // LBMDIM==3
+#              endif // LBMDIM==3
 
                if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) {
                        // normal ok and usable...
@@ -712,16 +762,27 @@ LbmFsgrSolver::mainLoop(int lev)
                // calculate macroscopic cell values
                LbmFloat oldUx, oldUy, oldUz;
                LbmFloat oldRho; // OLD rho = ccel->rho;
-#if OPT3D==0
-                       oldRho=RAC(ccel,0);
-                       oldUx = oldUy = oldUz = 0.0;
-                       for(int l=1; l<this->cDfNum; l++) {
-                               oldRho += RAC(ccel,l);
-                               oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
-                               oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
-                               oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
-                       } 
-#else // OPT3D==0
+#              define REFERENCE_PRESSURE 1.0 // always atmosphere...
+#              if OPT3D==0
+               oldRho=RAC(ccel,0);
+               oldUx = oldUy = oldUz = 0.0;
+               for(int l=1; l<this->cDfNum; l++) {
+                       oldRho += RAC(ccel,l);
+                       oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
+                       oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
+                       oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
+               } 
+               // reconstruct dist funcs from empty cells
+               FORDF1 {
+                       if(recons[ l ]) {
+                               m[ this->dfInv[l] ] = 
+                                       this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
+                                       this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
+                                       - MYDF( l );
+                       }
+               }
+               usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
+#              else // OPT3D==0
                oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
                                + RAC(ccel,dS ) + RAC(ccel,dE )
                                + RAC(ccel,dW ) + RAC(ccel,dT )
@@ -733,39 +794,25 @@ LbmFsgrSolver::mainLoop(int lev)
                                + RAC(ccel,dEB) + RAC(ccel,dWT)
                                + RAC(ccel,dWB);
 
-                       oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
+               oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
                                + RAC(ccel,dNE) - RAC(ccel,dNW)
                                + RAC(ccel,dSE) - RAC(ccel,dSW)
                                + RAC(ccel,dET) + RAC(ccel,dEB)
                                - RAC(ccel,dWT) - RAC(ccel,dWB);
 
-                       oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
+               oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
                                + RAC(ccel,dNE) + RAC(ccel,dNW)
                                - RAC(ccel,dSE) - RAC(ccel,dSW)
                                + RAC(ccel,dNT) + RAC(ccel,dNB)
                                - RAC(ccel,dST) - RAC(ccel,dSB);
 
-                       oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
+               oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
                                + RAC(ccel,dNT) - RAC(ccel,dNB)
                                + RAC(ccel,dST) - RAC(ccel,dSB)
                                + RAC(ccel,dET) - RAC(ccel,dEB)
                                + RAC(ccel,dWT) - RAC(ccel,dWB);
-#endif
 
                // now reconstruction
-#define REFERENCE_PRESSURE 1.0 // always atmosphere...
-#if OPT3D==0
-               // NOW - construct dist funcs from empty cells
-               FORDF1 {
-                       if(recons[ l ]) {
-                               m[ this->dfInv[l] ] = 
-                                       this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
-                                       this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
-                                       - MYDF( l );
-                       }
-               }
-               usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
-#else
                ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
                rho = REFERENCE_PRESSURE;
                usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
@@ -787,7 +834,7 @@ LbmFsgrSolver::mainLoop(int lev)
                if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
                if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
                if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
-#endif         
+#              endif           
 
 
                // inflow bc handling
@@ -835,9 +882,6 @@ LbmFsgrSolver::mainLoop(int lev)
                                && (this->mPartGenProb>0.0)) {
                        bool doAdd = true;
                        bool bndOk=true;
-                       //if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)||
-                                       //(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)||
-                                       //(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; }
                        if( (i<cutMin)||(i>this->mSizex-cutMin)||
                                        (j<cutMin)||(j>this->mSizey-cutMin)||
                                        (k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
@@ -853,27 +897,14 @@ LbmFsgrSolver::mainLoop(int lev)
                        LbmFloat prob = (rand()/(RAND_MAX+1.0));
                        LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl;
 
-                       // reduce probability in outer region
-                       //if( (i<cutConst)||(i>this->mSizex-cutConst) ){ prob *= 0.5; }
-                       //if( (j<cutConst)||(j>this->mSizey-cutConst) ){ prob *= 0.5; }
-                       //if( (k<cutConst)||(k>this->mSizez-cutConst) ){ prob *= 0.5; }
-                       // NEW TEST 040627
-                       //if( (i<cutConst)||(i>this->mSizex-cutConst) ){ doAdd=false; }
-                       //if( (j<cutConst)||(j>this->mSizey-cutConst) ){ doAdd=false; }
-                       //if( (k<cutConst)||(k>this->mSizez-cutConst) ){ doAdd=false; }
+                       // reduce probability in outer region?
                        const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
                        const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
                        LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
                        LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
                        if(pifac<0.) pifac=0.;
                        if(pjfac<0.) pjfac=0.;
-                       //const LbmFloat pkfac = 1.-(LbmFloat)(ABS(k-mLevel[mMaxRefine].lSizez/2))/(LbmFloat)(mLevel[mMaxRefine].lSizez/2);
-                       //errMsg("PROBTTT"," at "<<PRINT_IJK<<" prob"<<prob<<" pifac"<<pifac<<" pjfac"<<pjfac<<" "<<(basethresh*rl*pifac*pjfac));
-                       //prob *= pifac*pjfac; //*pkfac;
 
-//#define RWVEL_THRESH 1.0
-#define RWVEL_THRESH 1.5
-#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
                        //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
                        if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
                                // add
@@ -881,24 +912,12 @@ LbmFsgrSolver::mainLoop(int lev)
                                doAdd = false; // dont...
                        }
 
-//#define SLOWDOWNREGION (2*mCutoff)
-#if LBMDIM==3
-                       // normal
-#define SLOWDOWNREGION (this->mSizez/4)
-#else // LBMDIM==2
-                       // off
-#define SLOWDOWNREGION 10 
-#endif // LBMDIM==2
-#define P_LCSMQO 0.01
-
                        // "wind" disturbance
                        // use realworld relative velocity here instead?
                        if( (doAdd && 
                                        ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
                                        ||(k>this->mSizez-SLOWDOWNREGION)   ) {
                                LbmFloat nuz = uz;
-                               //? LbmFloat jdf; // = 0.05 * (rand()/(RAND_MAX+1.0));
-                               //? if(rl>RWVEL_WINDTHRESH) jdf *= (rl-RWVEL_WINDTHRESH);
                                if(k>this->mSizez-SLOWDOWNREGION) {
                                        // special case
                                        LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) );
@@ -917,7 +936,6 @@ LbmFsgrSolver::mainLoop(int lev)
                                        const LbmFloat add =  this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
                                        RAC(tcel,l) += add; }
                                }
-                               //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<PRINT_VEC(ux,uy,nuz)<<" rwv"<<PRINT_VEC(rux,ruy,ruz) );
                                //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
                        } // wind disturbance
 
@@ -972,17 +990,13 @@ LbmFsgrSolver::mainLoop(int lev)
                                mpParticles->getLast()->setType(type);
                                //if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
                                mpParticles->getLast()->setSize(size);
-                               //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<PRINT_VEC(ux,uy,uz) <<" RWu="<<PRINT_VEC(rux,ruy,ruz)<<" add"<<doAdd<<" pvel="<<pv );
                                //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
-                               //if(!bndOk){ mass -= val2fac*size*0.02; } // FORCEDISSOLVE
-                               //const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor!
-                               //mass -= val2fac*size*0.0015; // NTEST!
-                               //mass -= val2fac*size*0.001; // NTEST!
+                               //mass -= size*0.001; // NTEST!
                                mass -= size*0.0020; // NTEST!
-#if LBMDIM==2
+#                              if LBMDIM==2
                                mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
                                mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
-#endif // LBMDIM==2
+#                              endif // LBMDIM==2
                //errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel()  );
                                } // multiple parts
                        } // doAdd
@@ -1004,8 +1018,8 @@ LbmFsgrSolver::mainLoop(int lev)
                }
 
                // looks much nicer... LISTTRICK
-#if FSGR_LISTTRICK==1
-               if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
+#              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...
@@ -1023,7 +1037,7 @@ LbmFsgrSolver::mainLoop(int lev)
                                } 
                        }
                } // nobndfluid only */
-#endif
+#              endif
                //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
                
 
@@ -1032,11 +1046,6 @@ LbmFsgrSolver::mainLoop(int lev)
                        LbmPoint filledp; filledp.flag=0;
                        if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
                        filledp.x = i; filledp.y = j; filledp.z = k;
-#if PARALLEL==1
-                       //calcListFull[id].push_back( filledp );
-#else // PARALLEL==1
-                       //mListFull.push_back( filledp );
-#endif // PARALLEL==1
                        LIST_FULL(filledp);
                        //this->mNumFilledCells++; // DEBUG
                        calcCellsFilled++;
@@ -1045,11 +1054,6 @@ LbmFsgrSolver::mainLoop(int lev)
                        LbmPoint emptyp; emptyp.flag=0;
                        if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
                        emptyp.x = i; emptyp.y = j; emptyp.z = k;
-#if PARALLEL==1
-                       //calcListEmpty[id].push_back( emptyp );
-#else // PARALLEL==1
-                       //mListEmpty.push_back( emptyp );
-#endif // PARALLEL==1
                        LIST_EMPTY(emptyp);
                        //this->mNumEmptiedCells++; // DEBUG
                        calcCellsEmptied++;
@@ -1087,38 +1091,133 @@ LbmFsgrSolver::mainLoop(int lev)
                calcCurrentVolume += RAC(tcel,dFfrac);
 
                // interface cell handling done...
-       } // i
-       int i=0; //dummy
-       ADVANCE_POINTERS(2);
-       } // j
 
-#if COMPRESSGRIDS==1
-#if PARALLEL==1
-       //frintf(stderr," (id=%d k=%d) ",id,k);
-# pragma omp barrier
+//#if PARALLEL==1
+//#include "paraloopend.h"
+       //GRID_REGION_END();
+//#else // PARALLEL==1
+       //GRID_LOOPREG_END();
+       //GRID_REGION_END();
+//#endif // PARALLEL==1
+#if PARALLEL!=1
+       GRID_LOOPREG_END();
+#else // PARALLEL==1
+#include "paraloopend.h" // = GRID_LOOPREG_END();
 #endif // PARALLEL==1
-#else // COMPRESSGRIDS==1
-       int i=0; //dummy
-       ADVANCE_POINTERS(mLevel[lev].lSizex*2);
-#endif // COMPRESSGRIDS==1
-  } // all cell loop k,j,i
-
-       } // main loop region
 
-       // write vars from parallel computations to class
-       //errMsg("DFINI"," maxr l"<<mMaxRefine<<" cm="<<calcCurrentMass<<" cv="<<calcCurrentVolume );
+       // write vars from computations to class
        mLevel[lev].lmass    = calcCurrentMass;
        mLevel[lev].lvolume  = calcCurrentVolume;
        this->mNumFilledCells  = calcCellsFilled;
        this->mNumEmptiedCells = calcCellsEmptied;
        this->mNumUsedCells = calcNumUsedCells;
+}
+
+
+
+void 
+LbmFsgrSolver::preinitGrids()
+{
+       const int lev = mMaxRefine;
+       const bool doReduce = false;
+       const int gridLoopBound=0;
+
+       // touch both grids
+       for(int s=0; s<2; s++) {
+       
+       GRID_REGION_INIT();
 #if PARALLEL==1
-       PARA_FINISH();
+#include "paraloopstart.h"
 #endif // PARALLEL==1
+       GRID_REGION_START();
+       GRID_LOOP_START();
+               FORDF0{ RAC(ccel,l) = 0.; }
+               *pFlagSrc =0;
+               *pFlagDst =0;
+               //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
+#if PARALLEL!=1
+       GRID_LOOPREG_END();
+#else // PARALLEL==1
+#include "paraloopend.h" // = GRID_LOOPREG_END();
+#endif // PARALLEL==1
+       //GRID_REGION_END();
+       // TEST! */
 
-       // check other vars...?
+       /* dummy remove warnings */ 
+       calcCurrentMass = calcCurrentVolume = 0.;
+       calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
+       
+       // change grid
+  mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
+  mLevel[mMaxRefine].setCurr   ^= 1;
+       }
 }
 
+void 
+LbmFsgrSolver::standingFluidPreinit()
+{
+       const int lev = mMaxRefine;
+       const bool doReduce = false;
+       const int gridLoopBound=1;
+
+       GRID_REGION_INIT();
+#if PARALLEL==1
+#include "paraloopstart.h"
+#endif // PARALLEL==1
+       GRID_REGION_START();
+
+       LbmFloat rho, ux,uy,uz, usqr; 
+       CellFlagType nbflag[LBM_DFNUM];
+       LbmFloat m[LBM_DFNUM];
+       LbmFloat lcsmqo;
+#      if OPT3D==1 
+       LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
+       CellFlagType nbored=0;
+#      endif // OPT3D==true 
+
+       GRID_LOOP_START();
+               //FORDF0{ RAC(ccel,l) = 0.; }
+               //*pFlagSrc =0;
+               //*pFlagDst =0;
+               //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
+               const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
+               if( (currFlag & (CFEmpty|CFBnd)) ) continue;
+               //ccel = RACPNT(lev, i,j,k,workSet); 
+               //tcel = RACPNT(lev, i,j,k,otherSet);
+
+               if( (currFlag & (CFInter)) ) {
+                       // copy all values
+                       for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
+                       continue;
+               }
+
+               if( (currFlag & CFNoBndFluid)) {
+                       OPTIMIZED_STREAMCOLLIDE;
+               } else {
+                       FORDF1 {
+                               nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
+                       } 
+                       DEFAULT_STREAM;
+                       //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
+                       DEFAULT_COLLIDEG(mLevel[lev].gravity);
+               }
+               for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
+#if PARALLEL!=1
+       GRID_LOOPREG_END();
+#else // PARALLEL==1
+#include "paraloopend.h" // = GRID_LOOPREG_END();
+#endif // PARALLEL==1
+       //GRID_REGION_END();
+       // TEST! */
+
+       /* dummy remove warnings */ 
+       calcCurrentMass = calcCurrentVolume = 0.;
+       calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
+       
+       // change grid
+  mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
+  mLevel[mMaxRefine].setCurr   ^= 1;
+}
 
 
 /******************************************************************************
index 1a1e2c011726ae25e0fde68b6a71b8d23b45ed4f..2c99a853c6240e5a6d8ab9f7812841b01de38876 100644 (file)
 #define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
 #endif // FSGR_STRICT_DEBUG==1
 
+#if LBM_INCLUDE_TESTSOLVERS!=1
 
+#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+       ux += (grav)[0]; \
+       uy += (grav)[1]; \
+       uz += (grav)[2]; 
 
+#define TEST_IF_CHECK 
 
+#else // LBM_INCLUDE_TESTSOLVERS!=1
+// defined in test.h
+
+#define NEWDIRVELMOTEST 0
+#if NEWDIRVELMOTEST==1
 // off for non testing
+#undef PRECOLLIDE_MODS
 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
        ux += (grav)[0]; \
        uy += (grav)[1]; \
-       uz += (grav)[2]; \
-
-#define TEST_IF_CHECK 
+       uz += (grav)[2];  \
+  { \
+               int lev = mMaxRefine, nomb=0; \
+               LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
+               for(int l=1; l<this->cDfNum; l++) {  \
+                       if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
+                               if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
+                                       nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
+                                       nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
+                                       bcnt += 1.; \
+                               }       else { \
+                                       nomb++; \
+                               } \
+                       }  \
+               } \
+               if((bcnt>0.)&&(nomb==0)) { \
+                       ux = nux/bcnt; \
+                       uy = nuy/bcnt; \
+                       uz = nuz/bcnt;  \
+               } \
+       } 
+#else // NEWDIRVELMOTEST==1
+// off for non testing
+#endif // NEWDIRVELMOTEST==1
 
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
 
        
 /******************************************************************************
 // target set
 #define TSET(l) mLevel[(l)].setOther
 
+// handle mov. obj 
+#if FSGR_STRICT_DEBUG==1
+
+#define  LBMDS_ADDMOV(linv,l)  \
+        \
+       if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
+        \
+       LbmFloat dte=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)-(mSimulationTime+this->mpParam->getTimestep()); \
+       if( ABS(dte)< 1e-15 ) { \
+       m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+       } else { \
+       const int sdx = i+this->dfVecX[linv], sdy = j+this->dfVecY[linv], sdz = k+this->dfVecZ[linv]; \
+        \
+       errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" "<<PRINT_VEC(sdx,sdy,sdz)<<" t="<<(mSimulationTime+this->mpParam->getTimestep())<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)<<" dte="<<dte); \
+       debugMarkCell(lev,sdx,sdy,sdz); \
+       } \
+       } \
+
+
+
+#else // FSGR_STRICT_DEBUG==1
+
+#define  LBMDS_ADDMOV(linv,l)  \
+        \
+        \
+       if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
+        \
+       m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+       } \
+
+
+
+#endif // !FSGR_STRICT_DEBUG==1
+
 // treatment of freeslip reflection
 // used both for OPT and nonOPT
-#define DEFAULT_STREAM_FREESLIP(l,invl,mnbf) \
-               /*const int inv_l = this->dfInv[l];*/ \
-               int nb1 = 0, nb2 = 0;   /* is neighbor in this direction an obstacle? */\
-               LbmFloat newval = 0.0; /* new value for m[l], differs for free/part slip */\
-               const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \
-               if(dz==0) { \
-                       nb1 = !(RFLAG(lev, i,   j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
-                       nb2 = !(RFLAG(lev, i+dx,j,   k, SRCS(lev))&(CFFluid|CFInter)); \
-                       if((nb1)&&(!nb2)) { \
-                               /* x reflection */\
-                               newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
-                       } else  \
-                       if((!nb1)&&(nb2)) { \
-                               /* y reflection */\
-                               newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
-                       } else { \
-                               /* normal no slip in all other cases */\
-                               newval = QCELL(lev, i,j,k,SRCS(lev), invl); \
-                       } \
-               } else /* z=0 */\
-               if(dy==0) { \
-                       nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
-                       nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \
-                       if((nb1)&&(!nb2)) { \
-                               /* x reflection */\
-                               newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
-                       } else  \
-                       if((!nb1)&&(nb2)) { \
-                               /* z reflection */\
-                               newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
-                       } else { \
-                               /* normal no slip in all other cases */\
-                               newval = ( QCELL(lev, i,j,k,SRCS(lev), invl) ); \
-                       } \
-                       /* end y=0 */ \
-               } else { \
-                       /* x=0 */\
-                       nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
-                       nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
-                       if((nb1)&&(!nb2)) { \
-                               /* y reflection */\
-                               newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
-                       } else  \
-                       if((!nb1)&&(nb2)) { \
-                               /* z reflection */\
-                               newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
-                       } else { \
-                               /* normal no slip in all other cases */\
-                               newval = ( QCELL(lev, i,j,k,SRCS(lev), invl) ); \
-                       } \
-               } \
-               if(mnbf & CFBndPartslip) { /* part slip interpolation */ \
-                       const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \
-                       m[l] = RAC(ccel, this->dfInv[l] ) * partv + newval * (1.0-partv); /* part slip */ \
-               } else {\
-                       m[l] = newval; /* normal free slip*/\
-               }\
+#define  DEFAULT_STREAM_FREESLIP(l,invl,mnbf)  \
+        \
+       int nb1 = 0, nb2 = 0; \
+       LbmFloat newval = 0.0; \
+       const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \
+        \
+        \
+        \
+       const LbmFloat movadd = ( \
+       ((nbflag[invl]&CFBndMoving)&&(!(nbflag[l]&CFBnd))) ? \
+       (QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l)) : 0.); \
+        \
+       if(dz==0) { \
+       nb1 = !(RFLAG(lev, i,   j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
+       nb2 = !(RFLAG(lev, i+dx,j,   k, SRCS(lev))&(CFFluid|CFInter)); \
+       if((nb1)&&(!nb2)) { \
+        \
+       newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
+       } else \
+       if((!nb1)&&(nb2)) { \
+        \
+       newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
+       } else { \
+        \
+       newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \
+       } \
+       } else \
+       if(dy==0) { \
+       nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
+       nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \
+       if((nb1)&&(!nb2)) { \
+        \
+       newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
+       } else \
+       if((!nb1)&&(nb2)) { \
+        \
+       newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
+       } else { \
+        \
+       newval = RAC(ccel, this->dfInv[l] )  +movadd /* */; \
+       } \
+        \
+       } else \
+        \
+       { \
+        \
+       nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
+       nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
+       if((nb1)&&(!nb2)) { \
+        \
+       newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
+       } else \
+       if((!nb1)&&(nb2)) { \
+        \
+       newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
+       } else { \
+        \
+       newval = RAC(ccel, this->dfInv[l] )  +movadd /* */; \
+       } \
+       } \
+        \
+       if(mnbf & CFBndPartslip) { \
+       const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \
+        \
+       m[l] = (RAC(ccel, this->dfInv[l] )  +movadd /* d *(1./1.) */ ) * partv + newval * (1.0-partv); \
+       } else { \
+       m[l] = newval; \
+       } \
+        \
+
+
+
 
 // complete default stream&collide, 2d/3d
 /* read distribution funtions of adjacent cells = sweep step */ 
 #define MARKCELLCHECK \
        debugMarkCell(lev,i,j,k); CAUSE_PANIC;
 #define STREAMCHECK(id,ni,nj,nk,nl) \
-       if((m[nl] < -1.0) || (m[nl]>1.0)) {\
+       if((!(m[nl] > -1.0) && (m[nl]<1.0)) ) {\
                errMsg("STREAMCHECK","ID"<<id<<" Invalid streamed DF nl"<<nl<<" value:"<<m[nl]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\
                                " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther)  ); \
                /*FORDF0{ errMsg("STREAMCHECK"," at "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); } */ \
                MARKCELLCHECK; \
+               m[nl] = dfEquil[nl]; /* REPAIR */ \
        }
 #define COLLCHECK \
        if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\
                errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \
                /*FORDF0{ errMsg("COLLCHECK"," at? "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); }*/ \
+               rho=ux=uy=uz= 0.; /* REPAIR */ \
                MARKCELLCHECK; \
        }
 #else
 #endif
 
 // careful ux,uy,uz need to be inited before!
+#define  DEFAULT_STREAM  \
+       m[dC] = RAC(ccel,dC); \
+       STREAMCHECK(1, i,j,k, dC); \
+       FORDF1 { \
+       CellFlagType nbf = nbflag[ this->dfInv[l] ]; \
+       if(nbf & CFBnd) { \
+       if(nbf & CFBndNoslip) { \
+        \
+       m[l] = RAC(ccel, this->dfInv[l] ); \
+       LBMDS_ADDMOV(this->dfInv[l],l); \
+       STREAMCHECK(2, i,j,k, l); \
+       } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \
+        \
+       if(l<=LBMDIM*2) { \
+       m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); \
+       LBMDS_ADDMOV(this->dfInv[l],l); \
+       } else { \
+       const int inv_l = this->dfInv[l]; \
+       DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \
+       } \
+        \
+       } \
+       else { \
+       errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \
+       } \
+       } else { \
+       m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+       STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
+       } \
+       } \
+
 
-#define DEFAULT_STREAM \
-               m[dC] = RAC(ccel,dC); \
-               STREAMCHECK(1, i,j,k, dC); \
-               FORDF1 { \
-                       CellFlagType nbf = nbflag[ this->dfInv[l] ];\
-                       if(nbf & CFBnd) { \
-                               if(nbf & CFBndNoslip) { \
-                                       /* no slip, default */ \
-                                       m[l] = RAC(ccel, this->dfInv[l] ); /* noslip */ \
-                                       STREAMCHECK(2, i,j,k, l); \
-                               } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \
-                                       /* free slip */ \
-                                       if(l<=LBMDIM*2) { \
-                                               m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); /* noslip for <dim*2 */ \
-                                       } else { \
-                                               const int inv_l = this->dfInv[l]; \
-                                               DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \
-                                       } /* l>2*dim free slip */ \
-                                       \
-                               } /* type reflect */\
-                               else {\
-                                       errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \
-                               } \
-                               if(nbf&CFBndMoving) m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); /* obs. speed*/ \
-                       } else { \
-                               m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
-                               STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
-                       } \
-               }   
 
 
 // careful ux,uy,uz need to be inited before!
-#define DEFAULT_COLLIDEG(grav) \
-                       this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \
-                       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
-                       FORDF0 { RAC(tcel,l) = m[l]; }   \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
-                       COLLCHECK;
-#define OPTIMIZED_STREAMCOLLIDE \
-                       m[0] = RAC(ccel,0); \
-                       FORDF1 { /* df0 is set later on... */ \
-                               /* FIXME CHECK INV ? */\
-                               if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC;  \
-                               } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
-                               STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
-                       }   \
-                       rho=m[0]; \
-                       /* ux = mLevel[lev].gravity[0]; uy = [1]; uz = mLevel[lev].gravity[2]; */ \
-                       this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo   ); \
-                       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
-                       FORDF0 { RAC(tcel,l) = m[l]; } \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
-                       COLLCHECK;
+#define  DEFAULT_COLLIDEG(grav)  \
+       this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \
+       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
+       FORDF0 { RAC(tcel,l) = m[l]; } \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       COLLCHECK; \
 
-#else  // 3D, opt OPT3D==true
 
 
-// handle mov. obj 
-#if FSGR_STRICT_DEBUG==1
-#define LBMDS_ADDMOV(linv,l) if(nbflag[linv]&CFBndMoving){ \
-               if(QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)== (mSimulationTime+this->mpParam->getTimestep()) ) { \
-                       m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); /* obs. speed*/ \
-               } else { \
-                       CAUSE_PANIC; errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" t="<<mSimulationTime<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)); \
-               } \
-       }
-#else // FSGR_STRICT_DEBUG==1
-#define LBMDS_ADDMOV(linv,l) if(nbflag[linv]&CFBndMoving){ m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); } /* obs. speed*/ 
-#endif // FSGR_STRICT_DEBUG==1
+#define  OPTIMIZED_STREAMCOLLIDE  \
+       m[0] = RAC(ccel,0); \
+       FORDF1 { \
+        \
+       if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC; \
+       } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
+       STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
+       } \
+       rho=m[0]; \
+       DEFAULT_COLLIDEG(mLevel[lev].gravity) \
+
+
+
+#define  OPTIMIZED_STREAMCOLLIDE___UNUSED  \
+        \
+       this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo   ); \
+       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
+       FORDF0 { RAC(tcel,l) = m[l]; } \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       COLLCHECK; \
+
+
+
+#else  // 3D, opt OPT3D==true
+
 
 // default stream opt3d add moving bc val
-#define DEFAULT_STREAM \
-               m[dC] = RAC(ccel,dC); \
-               /* explicit streaming */ \
-               if((!nbored & CFBnd)) { \
-                       /* no boundary near?, no real speed diff.? */ \
-                       m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
-                       m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
-                       m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
-                       m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
-                       m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
-                       m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
-               } else { \
-                       /* explicit streaming, normal velocity always zero for obstacles */ \
-                       if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \
-                       if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \
-                       if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \
-                       if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \
-                       if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \
-                       if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \
-                       \
-                       /* also treat free slip here */ \
-                       if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} LBMDS_ADDMOV(dSW,dNE); } else { m[dNE] = CSRC_NE; } \
-                       if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} LBMDS_ADDMOV(dSE,dNW); } else { m[dNW] = CSRC_NW; } \
-                       if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} LBMDS_ADDMOV(dNW,dSE); } else { m[dSE] = CSRC_SE; } \
-                       if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} LBMDS_ADDMOV(dNE,dSW); } else { m[dSW] = CSRC_SW; } \
-                       if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} LBMDS_ADDMOV(dSB,dNT); } else { m[dNT] = CSRC_NT; } \
-                       if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} LBMDS_ADDMOV(dST,dNB); } else { m[dNB] = CSRC_NB; } \
-                       if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} LBMDS_ADDMOV(dNB,dST); } else { m[dST] = CSRC_ST; } \
-                       if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} LBMDS_ADDMOV(dNT,dSB); } else { m[dSB] = CSRC_SB; } \
-                       if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} LBMDS_ADDMOV(dWB,dET); } else { m[dET] = CSRC_ET; } \
-                       if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} LBMDS_ADDMOV(dWT,dEB); } else { m[dEB] = CSRC_EB; } \
-                       if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} LBMDS_ADDMOV(dEB,dWT); } else { m[dWT] = CSRC_WT; } \
-                       if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} LBMDS_ADDMOV(dET,dWB); } else { m[dWB] = CSRC_WB; } \
-               } 
-
-
-
-#define COLL_CALCULATE_DFEQ(dstarray) \
-                       dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
-                       dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
-                       dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
-                       dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
-                       dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
-                       dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; 
-#define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \
-                       lcsmqadd  = (srcArray##NE - lcsmeq[ dNE ]); \
-                       lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
-                       lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
-                       lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
-                       lcsmqo = (lcsmqadd*    lcsmqadd); \
-                       lcsmqadd  = (srcArray##ET - lcsmeq[  dET ]); \
-                       lcsmqadd -= (srcArray##EB - lcsmeq[  dEB ]); \
-                       lcsmqadd -= (srcArray##WT - lcsmeq[  dWT ]); \
-                       lcsmqadd += (srcArray##WB - lcsmeq[  dWB ]); \
-                       lcsmqo += (lcsmqadd*    lcsmqadd); \
-                       lcsmqadd  = (srcArray##NT - lcsmeq[  dNT ]); \
-                       lcsmqadd -= (srcArray##NB - lcsmeq[  dNB ]); \
-                       lcsmqadd -= (srcArray##ST - lcsmeq[  dST ]); \
-                       lcsmqadd += (srcArray##SB - lcsmeq[  dSB ]); \
-                       lcsmqo += (lcsmqadd*    lcsmqadd); \
-                       lcsmqo *= 2.0; \
-                       lcsmqadd  = (srcArray##E  -  lcsmeq[ dE  ]); \
-                       lcsmqadd += (srcArray##W  -  lcsmeq[ dW  ]); \
-                       lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
-                       lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
-                       lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
-                       lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
-                       lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
-                       lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
-                       lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
-                       lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
-                       lcsmqo += (lcsmqadd*    lcsmqadd); \
-                       lcsmqadd  = (srcArray##N  -  lcsmeq[ dN  ]); \
-                       lcsmqadd += (srcArray##S  -  lcsmeq[ dS  ]); \
-                       lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
-                       lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
-                       lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
-                       lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
-                       lcsmqadd += (srcArray##NT  - lcsmeq[ dNT ]); \
-                       lcsmqadd += (srcArray##NB  - lcsmeq[ dNB ]); \
-                       lcsmqadd += (srcArray##ST  - lcsmeq[ dST ]); \
-                       lcsmqadd += (srcArray##SB  - lcsmeq[ dSB ]); \
-                       lcsmqo += (lcsmqadd*    lcsmqadd); \
-                       lcsmqadd  = (srcArray##T  -  lcsmeq[ dT  ]); \
-                       lcsmqadd += (srcArray##B  -  lcsmeq[ dB  ]); \
-                       lcsmqadd += (srcArray##NT -  lcsmeq[ dNT ]); \
-                       lcsmqadd += (srcArray##NB -  lcsmeq[ dNB ]); \
-                       lcsmqadd += (srcArray##ST -  lcsmeq[ dST ]); \
-                       lcsmqadd += (srcArray##SB -  lcsmeq[ dSB ]); \
-                       lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
-                       lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
-                       lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
-                       lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
-                       lcsmqo += (lcsmqadd*    lcsmqadd); \
-                       lcsmqo = sqrt(lcsmqo); /* FIXME check effect of sqrt*/ \
+#define  DEFAULT_STREAM  \
+       m[dC] = RAC(ccel,dC); \
+        \
+       if((!nbored & CFBnd)) { \
+        \
+       m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
+       m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
+       m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
+       m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
+       m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
+       m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
+       } else { \
+        \
+       if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \
+       if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \
+       if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \
+       if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \
+       if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \
+       if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \
+        \
+        \
+       if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} } else { m[dNE] = CSRC_NE; } \
+       if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} } else { m[dNW] = CSRC_NW; } \
+       if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} } else { m[dSE] = CSRC_SE; } \
+       if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} } else { m[dSW] = CSRC_SW; } \
+       if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} } else { m[dNT] = CSRC_NT; } \
+       if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} } else { m[dNB] = CSRC_NB; } \
+       if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} } else { m[dST] = CSRC_ST; } \
+       if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} } else { m[dSB] = CSRC_SB; } \
+       if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} } else { m[dET] = CSRC_ET; } \
+       if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} } else { m[dEB] = CSRC_EB; } \
+       if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} } else { m[dWT] = CSRC_WT; } \
+       if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} } else { m[dWB] = CSRC_WB; } \
+       } \
+
+
+
+
+
+#define  COLL_CALCULATE_DFEQ(dstarray)  \
+       dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
+       dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
+       dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
+       dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
+       dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
+       dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; \
+
+
+
+#define  COLL_CALCULATE_NONEQTENSOR(csolev, srcArray )  \
+       lcsmqadd  = (srcArray##NE - lcsmeq[ dNE ]); \
+       lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
+       lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
+       lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
+       lcsmqo = (lcsmqadd*    lcsmqadd); \
+       lcsmqadd  = (srcArray##ET - lcsmeq[  dET ]); \
+       lcsmqadd -= (srcArray##EB - lcsmeq[  dEB ]); \
+       lcsmqadd -= (srcArray##WT - lcsmeq[  dWT ]); \
+       lcsmqadd += (srcArray##WB - lcsmeq[  dWB ]); \
+       lcsmqo += (lcsmqadd*    lcsmqadd); \
+       lcsmqadd  = (srcArray##NT - lcsmeq[  dNT ]); \
+       lcsmqadd -= (srcArray##NB - lcsmeq[  dNB ]); \
+       lcsmqadd -= (srcArray##ST - lcsmeq[  dST ]); \
+       lcsmqadd += (srcArray##SB - lcsmeq[  dSB ]); \
+       lcsmqo += (lcsmqadd*    lcsmqadd); \
+       lcsmqo *= 2.0; \
+       lcsmqadd  = (srcArray##E  -  lcsmeq[ dE  ]); \
+       lcsmqadd += (srcArray##W  -  lcsmeq[ dW  ]); \
+       lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
+       lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
+       lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
+       lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
+       lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
+       lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
+       lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
+       lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
+       lcsmqo += (lcsmqadd*    lcsmqadd); \
+       lcsmqadd  = (srcArray##N  -  lcsmeq[ dN  ]); \
+       lcsmqadd += (srcArray##S  -  lcsmeq[ dS  ]); \
+       lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
+       lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
+       lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
+       lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
+       lcsmqadd += (srcArray##NT  - lcsmeq[ dNT ]); \
+       lcsmqadd += (srcArray##NB  - lcsmeq[ dNB ]); \
+       lcsmqadd += (srcArray##ST  - lcsmeq[ dST ]); \
+       lcsmqadd += (srcArray##SB  - lcsmeq[ dSB ]); \
+       lcsmqo += (lcsmqadd*    lcsmqadd); \
+       lcsmqadd  = (srcArray##T  -  lcsmeq[ dT  ]); \
+       lcsmqadd += (srcArray##B  -  lcsmeq[ dB  ]); \
+       lcsmqadd += (srcArray##NT -  lcsmeq[ dNT ]); \
+       lcsmqadd += (srcArray##NB -  lcsmeq[ dNB ]); \
+       lcsmqadd += (srcArray##ST -  lcsmeq[ dST ]); \
+       lcsmqadd += (srcArray##SB -  lcsmeq[ dSB ]); \
+       lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
+       lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
+       lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
+       lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
+       lcsmqo += (lcsmqadd*    lcsmqadd); \
+       lcsmqo = sqrt(lcsmqo); \
+
+
 
 //                     COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega); 
 
 // careful - need lcsmqo 
-#define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \
-                       dstomega =  1.0/\
-                                       ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*(\
-                                                       -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
-                                                       / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
-                                               ) +0.5 ); 
-
-#define DEFAULT_COLLIDE_LES(grav) \
-                       rho = + MSRC_C  + MSRC_N  \
-                               + MSRC_S  + MSRC_E  \
-                               + MSRC_W  + MSRC_T  \
-                               + MSRC_B  + MSRC_NE \
-                               + MSRC_NW + MSRC_SE \
-                               + MSRC_SW + MSRC_NT \
-                               + MSRC_NB + MSRC_ST \
-                               + MSRC_SB + MSRC_ET \
-                               + MSRC_EB + MSRC_WT \
-                               + MSRC_WB; \
-                       \
-                       ux = MSRC_E - MSRC_W \
-                               + MSRC_NE - MSRC_NW \
-                               + MSRC_SE - MSRC_SW \
-                               + MSRC_ET + MSRC_EB \
-                               - MSRC_WT - MSRC_WB ;  \
-                       \
-                       uy = MSRC_N - MSRC_S \
-                               + MSRC_NE + MSRC_NW \
-                               - MSRC_SE - MSRC_SW \
-                               + MSRC_NT + MSRC_NB \
-                               - MSRC_ST - MSRC_SB ;  \
-                       \
-                       uz = MSRC_T - MSRC_B \
-                               + MSRC_NT - MSRC_NB \
-                               + MSRC_ST - MSRC_SB \
-                               + MSRC_ET - MSRC_EB \
-                               + MSRC_WT - MSRC_WB ;  \
-                       PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
-                       COLL_CALCULATE_DFEQ(lcsmeq); \
-                       COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \
-                       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
-                       CSMOMEGA_STATS(lev,lcsmomega); \
-                       \
-                       RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
-                       \
-                       RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
-                       RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
-                       RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
-                       RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
-                       RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
-                       RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
-                       \
-                       RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
-                       RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
-                       RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
-                       RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
-                       RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
-                       RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
-                       RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
-                       RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
-                       RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
-                       RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
-                       RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
-                       RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; 
-
-#define DEFAULT_COLLIDE_NOLES(grav) \
-                       rho = + MSRC_C  + MSRC_N  \
-                               + MSRC_S  + MSRC_E  \
-                               + MSRC_W  + MSRC_T  \
-                               + MSRC_B  + MSRC_NE \
-                               + MSRC_NW + MSRC_SE \
-                               + MSRC_SW + MSRC_NT \
-                               + MSRC_NB + MSRC_ST \
-                               + MSRC_SB + MSRC_ET \
-                               + MSRC_EB + MSRC_WT \
-                               + MSRC_WB; \
-                       \
-                       ux = MSRC_E - MSRC_W \
-                               + MSRC_NE - MSRC_NW \
-                               + MSRC_SE - MSRC_SW \
-                               + MSRC_ET + MSRC_EB \
-                               - MSRC_WT - MSRC_WB ;  \
-                       \
-                       uy = MSRC_N - MSRC_S \
-                               + MSRC_NE + MSRC_NW \
-                               - MSRC_SE - MSRC_SW \
-                               + MSRC_NT + MSRC_NB \
-                               - MSRC_ST - MSRC_SB ;  \
-                       \
-                       uz = MSRC_T - MSRC_B \
-                               + MSRC_NT - MSRC_NB \
-                               + MSRC_ST - MSRC_SB \
-                               + MSRC_ET - MSRC_EB \
-                               + MSRC_WT - MSRC_WB ;  \
-                       PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
-                       \
-                       RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C  + OMEGA(lev)*EQC ; \
-                       \
-                       RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N  + OMEGA(lev)*EQN ; \
-                       RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S  + OMEGA(lev)*EQS ; \
-                       RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E  + OMEGA(lev)*EQE ; \
-                       RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W  + OMEGA(lev)*EQW ; \
-                       RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T  + OMEGA(lev)*EQT ; \
-                       RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B  + OMEGA(lev)*EQB ; \
-                       \
-                       RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
-                       RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
-                       RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
-                       RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
-                       RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
-                       RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
-                       RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
-                       RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
-                       RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
-                       RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
-                       RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
-                       RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; 
-
-
-
-#define OPTIMIZED_STREAMCOLLIDE_LES \
-                       /* only surrounded by fluid cells...!, so safe streaming here... */ \
-                       m[dC ] = CSRC_C ; \
-                       m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
-                       m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
-                       m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
-                       m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
-                       m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
-                       m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
-                       \
-                       rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  \
-                               + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
-                               + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
-                       ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
-                               + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB;  \
-                       uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
-                               + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB;  \
-                       uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
-                               + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB;  \
-                       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
-                       COLL_CALCULATE_DFEQ(lcsmeq); \
-                       COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
-                       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
-                       CSMOMEGA_STATS(lev,lcsmomega); \
-                       \
-                       RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
-                       RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ];  \
-                       RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ];  \
-                       RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
-                       RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ];  \
-                       RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ];  \
-                       RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
-                       \
-                       RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE];  \
-                       RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW];  \
-                       RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE];  \
-                       RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
-                       \
-                       RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT];  \
-                       RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB];  \
-                       RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST];  \
-                       RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
-                       \
-                       RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET];  \
-                       RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
-                       RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT];  \
-                       RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB];  \
-
-#define OPTIMIZED_STREAMCOLLIDE_UNUSED \
-                       /* only surrounded by fluid cells...!, so safe streaming here... */ \
-                       rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
-                               + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
-                               + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
-                       ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
-                               + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB;  \
-                       uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
-                               + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB;  \
-                       uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
-                               + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB;  \
-                       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
-                       COLL_CALCULATE_DFEQ(lcsmeq); \
-                       COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
-                       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
-                       \
-                       RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C  + lcsmomega*EQC ; \
-                       RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N  + lcsmomega*lcsmeq[ dN ];  \
-                       RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S  + lcsmomega*lcsmeq[ dS ];  \
-                       RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E  + lcsmomega*lcsmeq[ dE ]; \
-                       RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W  + lcsmomega*lcsmeq[ dW ];  \
-                       RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T  + lcsmomega*lcsmeq[ dT ];  \
-                       RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B  + lcsmomega*lcsmeq[ dB ]; \
-                       \
-                       RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE];  \
-                       RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW];  \
-                       RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE];  \
-                       RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
-                       \
-                       RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT];  \
-                       RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB];  \
-                       RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST];  \
-                       RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
-                       \
-                       RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET];  \
-                       RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
-                       RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT];  \
-                       RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB];  \
-
-#define OPTIMIZED_STREAMCOLLIDE_NOLES \
-                       /* only surrounded by fluid cells...!, so safe streaming here... */ \
-                       rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
-                               + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
-                               + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
-                       ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
-                               + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB;  \
-                       uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
-                               + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB;  \
-                       uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
-                               + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB;  \
-                       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
-                       RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C  + OMEGA(lev)*EQC ; \
-                       RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N  + OMEGA(lev)*EQN ;  \
-                       RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S  + OMEGA(lev)*EQS ;  \
-                       RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E  + OMEGA(lev)*EQE ; \
-                       RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W  + OMEGA(lev)*EQW ;  \
-                       RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T  + OMEGA(lev)*EQT ;  \
-                       RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B  + OMEGA(lev)*EQB ; \
-                        \
-                       RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE;  \
-                       RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW;  \
-                       RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE;  \
-                       RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
-                        \
-                       RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT;  \
-                       RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB;  \
-                       RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST;  \
-                       RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
-                        \
-                       RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET;  \
-                       RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
-                       RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT;  \
-                       RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB;  \
-
-
-// debug version1
-#define STREAMCHECK(ni,nj,nk,nl) 
-#define COLLCHECK
-#define OPTIMIZED_STREAMCOLLIDE_DEBUG \
-                       m[0] = RAC(ccel,0); \
-                       FORDF1 { /* df0 is set later on... */ \
-                               if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC;\
-                               } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
-                               STREAMCHECK(9, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
-                       }   \
-                       /* rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; */ \
-                       this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo   ); \
-                       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
-                       FORDF0 { RAC(tcel,l) = m[l]; } \
-                       usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
-                       COLLCHECK;
-
-
-
-// more debugging
-/*DEBUG \
-                       m[0] = RAC(ccel,0); \
-                       FORDF1 { \
-                               if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl");CAUSE_PANIC; \
-                               } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
-                       }   \
-errMsg("T","QSDM at %d,%d,%d  lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lcsmomega ); \
-                       rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
-                       this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago  , &mDebugOmegaRet, &lcsmqo  ); \
-                       CSMOMEGA_STATS(lev,mDebugOmegaRet); \
-                       */
+#define  COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega )  \
+       dstomega =  1.0/ \
+       ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*( \
+       -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
+       / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
+       ) +0.5 ); \
+
+
+
+#define  DEFAULT_COLLIDE_LES(grav)  \
+       rho = + MSRC_C  + MSRC_N \
+       + MSRC_S  + MSRC_E \
+       + MSRC_W  + MSRC_T \
+       + MSRC_B  + MSRC_NE \
+       + MSRC_NW + MSRC_SE \
+       + MSRC_SW + MSRC_NT \
+       + MSRC_NB + MSRC_ST \
+       + MSRC_SB + MSRC_ET \
+       + MSRC_EB + MSRC_WT \
+       + MSRC_WB; \
+        \
+       ux = MSRC_E - MSRC_W \
+       + MSRC_NE - MSRC_NW \
+       + MSRC_SE - MSRC_SW \
+       + MSRC_ET + MSRC_EB \
+       - MSRC_WT - MSRC_WB ; \
+        \
+       uy = MSRC_N - MSRC_S \
+       + MSRC_NE + MSRC_NW \
+       - MSRC_SE - MSRC_SW \
+       + MSRC_NT + MSRC_NB \
+       - MSRC_ST - MSRC_SB ; \
+        \
+       uz = MSRC_T - MSRC_B \
+       + MSRC_NT - MSRC_NB \
+       + MSRC_ST - MSRC_SB \
+       + MSRC_ET - MSRC_EB \
+       + MSRC_WT - MSRC_WB ; \
+       PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       COLL_CALCULATE_DFEQ(lcsmeq); \
+       COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \
+       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+       CSMOMEGA_STATS(lev,lcsmomega); \
+        \
+       RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
+        \
+       RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
+       RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
+       RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
+       RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
+       RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
+       RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
+        \
+       RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+       RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+       RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+       RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+       RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+       RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+       RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
+       RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+       RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
+       RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+       RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+       RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define  DEFAULT_COLLIDE_NOLES(grav)  \
+       rho = + MSRC_C  + MSRC_N \
+       + MSRC_S  + MSRC_E \
+       + MSRC_W  + MSRC_T \
+       + MSRC_B  + MSRC_NE \
+       + MSRC_NW + MSRC_SE \
+       + MSRC_SW + MSRC_NT \
+       + MSRC_NB + MSRC_ST \
+       + MSRC_SB + MSRC_ET \
+       + MSRC_EB + MSRC_WT \
+       + MSRC_WB; \
+        \
+       ux = MSRC_E - MSRC_W \
+       + MSRC_NE - MSRC_NW \
+       + MSRC_SE - MSRC_SW \
+       + MSRC_ET + MSRC_EB \
+       - MSRC_WT - MSRC_WB ; \
+        \
+       uy = MSRC_N - MSRC_S \
+       + MSRC_NE + MSRC_NW \
+       - MSRC_SE - MSRC_SW \
+       + MSRC_NT + MSRC_NB \
+       - MSRC_ST - MSRC_SB ; \
+        \
+       uz = MSRC_T - MSRC_B \
+       + MSRC_NT - MSRC_NB \
+       + MSRC_ST - MSRC_SB \
+       + MSRC_ET - MSRC_EB \
+       + MSRC_WT - MSRC_WB ; \
+       PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+        \
+       RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C  + OMEGA(lev)*EQC ; \
+        \
+       RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N  + OMEGA(lev)*EQN ; \
+       RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S  + OMEGA(lev)*EQS ; \
+       RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E  + OMEGA(lev)*EQE ; \
+       RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W  + OMEGA(lev)*EQW ; \
+       RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T  + OMEGA(lev)*EQT ; \
+       RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B  + OMEGA(lev)*EQB ; \
+        \
+       RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
+       RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
+       RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
+       RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
+       RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
+       RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
+       RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
+       RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
+       RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
+       RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
+       RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
+       RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; \
+
+
+
+
+
+#define  OPTIMIZED_STREAMCOLLIDE_LES  \
+        \
+       m[dC ] = CSRC_C ; \
+       m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
+       m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
+       m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
+       m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
+       m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
+       m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
+        \
+       rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T \
+       + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
+       + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
+       ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
+       + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; \
+       uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
+       + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; \
+       uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
+       + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; \
+       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       COLL_CALCULATE_DFEQ(lcsmeq); \
+       COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
+       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+       CSMOMEGA_STATS(lev,lcsmomega); \
+        \
+       RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
+       RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
+       RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
+       RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
+       RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
+       RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
+       RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
+        \
+       RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+       RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+       RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+       RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+        \
+       RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+       RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+       RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
+       RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+        \
+       RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
+       RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+       RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+       RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define  OPTIMIZED_STREAMCOLLIDE_UNUSED  \
+        \
+       rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T \
+       + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
+       + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
+       ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
+       + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
+       uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
+       + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
+       uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
+       + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
+       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       COLL_CALCULATE_DFEQ(lcsmeq); \
+       COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
+       COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+        \
+       RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C  + lcsmomega*EQC ; \
+       RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N  + lcsmomega*lcsmeq[ dN ]; \
+       RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S  + lcsmomega*lcsmeq[ dS ]; \
+       RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E  + lcsmomega*lcsmeq[ dE ]; \
+       RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W  + lcsmomega*lcsmeq[ dW ]; \
+       RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T  + lcsmomega*lcsmeq[ dT ]; \
+       RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B  + lcsmomega*lcsmeq[ dB ]; \
+        \
+       RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+       RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+       RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+       RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+        \
+       RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+       RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+       RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST]; \
+       RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+        \
+       RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET]; \
+       RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+       RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+       RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define  OPTIMIZED_STREAMCOLLIDE_NOLES  \
+        \
+       rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T \
+       + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
+       + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
+       ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
+       + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
+       uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
+       + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
+       uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
+       + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
+       PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+       usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+       RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C  + OMEGA(lev)*EQC ; \
+       RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N  + OMEGA(lev)*EQN ; \
+       RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S  + OMEGA(lev)*EQS ; \
+       RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E  + OMEGA(lev)*EQE ; \
+       RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W  + OMEGA(lev)*EQW ; \
+       RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T  + OMEGA(lev)*EQT ; \
+       RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B  + OMEGA(lev)*EQB ; \
+        \
+       RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE; \
+       RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW; \
+       RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE; \
+       RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
+        \
+       RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT; \
+       RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB; \
+       RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST; \
+       RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
+        \
+       RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET; \
+       RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
+       RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT; \
+       RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB; \
+
+
+
+
+
+// LES switching for OPT3D
 #if USE_LES==1
 #define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_LES(grav)
 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES
@@ -723,7 +797,7 @@ errMsg("T","QSDM at %d,%d,%d  lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lc
 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES
 #endif
 
-#endif
+#endif  // 3D, opt OPT3D==true
 
 #define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz,  CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \
                        if(Cusqr>CmMaxVlen) { \
@@ -779,10 +853,6 @@ errMsg("T","QSDM at %d,%d,%d  lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lc
                                }
                                // end ADD_INT_DFSCHECK
                                
-                               //if(   !(RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) & CFUnused) ){
-                               //errMsg("INTFLAGUNU", PRINT_VEC(i,j,k)<<" child at "<<PRINT_VEC((ix),(iy),(iz)) );
-                               //if(iy==15) errMsg("IFFC", PRINT_VEC(i,j,k)<<" child interpolated at "<<PRINT_VEC((ix),(iy),(iz)) );
-                               //if(((ix)>10)&&(iy>5)&&(iz>5)) { debugMarkCell(lev+1, (ix),(iy),(iz) ); }
 #define INTUNUTCHECK(ix,iy,iz) \
                                if(     (RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) != (CFFluid|CFGrFromCoarse)) ){\
                                        errMsg("INTFLAGUNU_CHECK", PRINT_VEC(i,j,k)<<" child not unused at l"<<(lev+1)<<" "<<PRINT_VEC((ix),(iy),(iz))<<" flag: "<<  RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) ); \
@@ -1078,7 +1148,7 @@ inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmF
 
 // "normal" collision
 inline void LbmFsgrSolver::collideArrays(
-               int i, int j, int k, // position - more for debugging
+               int lev, int i, int j, int k, // position - more for debugging
                LbmFloat df[],                          
                LbmFloat &outrho, // out only!
                // velocity modifiers (returns actual velocity!)
@@ -1104,6 +1174,7 @@ inline void LbmFsgrSolver::collideArrays(
                uz  += (this->dfDvecZ[l]*df[l]);  
        }  
 
+
        PRECOLLIDE_MODS(rho,ux,uy,uz, gravity);
        for(l=0; l<this->cDfNum; l++) { 
                feq[l] = getCollideEq(l,rho,ux,uy,uz); 
index a38307a3b02a8eb8aa8e0936864fd17285574f61..6b7f6f2ed9c83bc536bc19024bb3b6ef4f62d4ae 100644 (file)
 #include "solver_relax.h"
 #include "particletracer.h"
 
+// MPT
+#include "ntl_world.h"
+#include "simulation_object.h"
+
 
 /******************************************************************************
  * helper functions
  *****************************************************************************/
 
+// try to enhance surface?
+#define SURFACE_ENH 2
+
 //! for raytracing
 void LbmFsgrSolver::prepareVisualization( void ) {
        int lev = mMaxRefine;
        int workSet = mLevel[lev].setCurr;
 
+       int mainGravDir=0;
+       LbmFloat mainGravLen = 0.;
+       FORDF1{
+               LbmFloat thisGravLen = dot(LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l]), getNormalized(mLevel[mMaxRefine].gravity) );
+               if(thisGravLen>mainGravLen) {
+                       mainGravLen = thisGravLen;
+                       mainGravDir = l;
+               }
+               //errMsg("GRAVT","l"<<l<<" dfv"<<LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l])<<" g"<< mLevel[mMaxRefine].gravity<< " mgl"<<mainGravLen<<" mgd"<<mainGravDir);
+       }
+
        //make same prepareVisualization and getIsoSurface...
 #if LBMDIM==2
        // 2d, place in the middle of isofield slice (k=2)
@@ -45,8 +63,17 @@ void LbmFsgrSolver::prepareVisualization( void ) {
        }
 #endif // LBMDIM==2
 
+       // MPT, ignore
+       //if( strstr( this->getName().c_str(), "mpfluid2" ) != NULL) {
+               //errMsg("MPT TESTX","skipping mpfluid2");
+               //return;
+       //}
+       if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
+               mpIso->resetAll(0.);
+       }
+
 
-       LbmFloat minval = this->mIsoValue*1.1; // / mIsoWeight[13]; 
+       LbmFloat minval = this->mIsoValue*1.05; // / mIsoWeight[13]; 
        // add up...
        float val = 0.0;
        for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 
@@ -54,60 +81,67 @@ void LbmFsgrSolver::prepareVisualization( void ) {
     for(int i=1;i<mLevel[lev].lSizex-1;i++) {
                        const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
                        //if(cflag&(CFBnd|CFEmpty)) {
-                       if(cflag&(CFBnd)) {
+
+#if SURFACE_ENH==0
+
+                       // no enhancements...
+                       if( (cflag&(CFFluid|CFUnused)) ) {
+                               val = 1.;
+                       } else if( (cflag&CFInter) ) {
+                               val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
+                       } else {
                                continue;
+                       }
 
-                       } else if( (cflag&CFEmpty) ) {
-                               //continue; // OFF DEBUG
+#else // SURFACE_ENH!=1
+                       if(cflag&CFBnd) {
+                               // treated in second loop
+                               continue;
+                       } else if(cflag&CFUnused) {
+                               val = 1.;
+                       } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid)) {
+                               // optimized fluid
+                               val = 1.;
+                       } else if( (cflag&(CFEmpty|CFInter|CFFluid)) ) {
                                int noslipbnd = 0;
                                int intercnt = 0;
-                               CellFlagType nbored;
                                FORDF1 { 
                                        const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
-                                       if((l==6)&&(nbflag&CFBnd)){ 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) ) {
+                                       if(l!=mainGravDir) continue; // only check bnd along main grav. dir
+                                       //if((nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
+                                       if((nbflag&CFBnd)){ noslipbnd=1; }
+                               }
 
-                       } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
-                               // optimized fluid
-                               val = 1.;
+                               if(cflag&CFEmpty) {
+                                       // special empty treatment
+                                       if((noslipbnd)&&(intercnt>6)) {
+                                               *this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
+                                       } else if((noslipbnd)&&(intercnt>0)) {
+                                               // necessary?
+                                               *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.9;
+                                       } else {
+                                               // nothing to do...
+                                       }
 
-                       } else if( (cflag&(CFInter|CFFluid)) ) {
-                               int noslipbnd = 0;
-                               FORDF1 { 
-                                       const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
-                                       if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; l=100; } 
-                               }
-                               // no empty nb interface cells are treated as full
-                               if(cflag&(CFNoNbEmpty|CFFluid)) {
+                                       continue;
+                               } else if(cflag&(CFNoNbEmpty|CFFluid)) {
+                                       // no empty nb interface cells are treated as full
                                        val=1.0;
+                               } else {
+                                       val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
                                }
-                               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] ); 
                                }
-                               // */
+#endif // SURFACE_ENH>0
 
-                       } else { // unused?
+                       } else { // all others, unused?
                                continue;
-                               // old fluid val = 1.0; 
-                       } // */
+                       } 
 
                        *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
                        *this->mpIso->lbmGetData( i   , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 
@@ -148,6 +182,105 @@ void LbmFsgrSolver::prepareVisualization( void ) {
                        *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
        }
 
+       // TEST!?
+#if SURFACE_ENH>=2
+
+       if(mFsSurfGenSetting&fssgNoObs) {
+               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++) {
+                               const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
+                               if(cflag&(CFBnd)) {
+                                       CellFlagType nbored=0;
+                                       LbmFloat avgfill=0.,avgfcnt=0.;
+                                       FORDF1 { 
+                                               const int ni = i+this->dfVecX[l];
+                                               const int nj = j+this->dfVecY[l];
+                                               const int nk = ZKOFF+this->dfVecZ[l];
+
+                                               const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
+                                               nbored |= nbflag;
+                                               if(nbflag&CFInter) {
+                                                       avgfill += QCELL(lev, ni,nj,nk, workSet,dFfrac); avgfcnt += 1.;
+                                               } else if(nbflag&CFFluid) {
+                                                       avgfill += 1.; avgfcnt += 1.;
+                                               } else if(nbflag&CFEmpty) {
+                                                       avgfill += 0.; avgfcnt += 1.;
+                                               }
+
+                                               //if( (ni<0) || (nj<0) || (nk<0) 
+                                                //|| (ni>=mLevel[mMaxRefine].lSizex) 
+                                                //|| (nj>=mLevel[mMaxRefine].lSizey) 
+                                                //|| (nk>=mLevel[mMaxRefine].lSizez) ) continue;
+                                       }
+
+                                       if(nbored&CFInter) {
+                                               if(avgfcnt>0.) avgfill/=avgfcnt;
+                                               *this->mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue;
+                                       } 
+                                       else if(nbored&CFFluid) {
+                                               *this->mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue;
+                                       }
+
+                               }
+                       }
+
+               // move surface towards inner "row" of obstacle
+               // cells if necessary (all obs cells without fluid/inter
+               // nbs (=iso==0) next to obstacles...)
+               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++) {
+                                       const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
+                                       if( (cflag&(CFBnd)) && (*this->mpIso->lbmGetData(i,j,ZKOFF)==0.)) {
+                                               int bndnbcnt=0;
+                                               FORDF1 { 
+                                                       const int ni = i+this->dfVecX[l];
+                                                       const int nj = j+this->dfVecY[l];
+                                                       const int nk = ZKOFF+this->dfVecZ[l];
+                                                       const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
+                                                       if(nbflag&CFBnd) bndnbcnt++;
+                                               }
+                                               if(bndnbcnt>0) *this->mpIso->lbmGetData(i,j,ZKOFF)=this->mIsoValue*0.95; 
+                                       }
+                               }
+       }
+       // */
+
+       if(mFsSurfGenSetting&fssgNoNorth) 
+               for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
+                       for(int j=0;j<mLevel[lev].lSizey-0;j++) {
+                               *this->mpIso->lbmGetData(0,                   j,ZKOFF) = *this->mpIso->lbmGetData(1,                   j,ZKOFF);
+                       }
+       if(mFsSurfGenSetting&fssgNoEast) 
+               for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
+                       for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+                               *this->mpIso->lbmGetData(i,0,                   ZKOFF) = *this->mpIso->lbmGetData(i,1,                   ZKOFF);
+                       }
+       if(mFsSurfGenSetting&fssgNoSouth) 
+               for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
+                       for(int j=0;j<mLevel[lev].lSizey-0;j++) {
+                               *this->mpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *this->mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF);
+                       }
+       if(mFsSurfGenSetting&fssgNoWest) 
+               for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
+                       for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+                               *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF);
+                       }
+       if(LBMDIM>2) {
+               if(mFsSurfGenSetting&fssgNoBottom) 
+                       for(int j=0;j<mLevel[lev].lSizey-0;j++) 
+                               for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+                                       *this->mpIso->lbmGetData(i,j,0                   ) = *this->mpIso->lbmGetData(i,j,1                   );
+                               } 
+               if(mFsSurfGenSetting&fssgNoTop) 
+                       for(int j=0;j<mLevel[lev].lSizey-0;j++) 
+                               for(int i=0;i<mLevel[lev].lSizex-0;i++) {
+                                       *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2);
+                               } 
+       }
+#endif // SURFACE_ENH>=2
+
 
        // update preview, remove 2d?
        if((this->mOutputSurfacePreview)&&(LBMDIM==3)) {
@@ -215,7 +348,47 @@ void LbmFsgrSolver::prepareVisualization( void ) {
                }
        }
 
-       // correction
+       // MPT
+#if LBM_INCLUDE_TESTSOLVERS==1
+       if( ( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) && (mInitDone) ) {
+               errMsg("MPT TESTX","special mpfluid1");
+               vector<SimulationObject*> *sims = mpGlob->getSims();
+               for(int simi=0; simi<(int)(sims->size()); simi++) {
+                       SimulationObject *sim = (*sims)[simi];
+                       if( strstr( sim->getName().c_str(), "mpfluid2" ) != NULL) {
+                               LbmFsgrSolver *fsgr = (LbmFsgrSolver *)(sim->getSolver());
+                               int msyLev = mMaxRefine;
+                               int simLev = fsgr->mMaxRefine;
+
+                               IsoSurface* simIso = fsgr->mpIso;
+                               int idst = mLevel[msyLev].lSizex-2 -1; // start at boundary
+                               int isrc = 0                    +2;
+                               if((simIso)&&(fsgr->mInitDone)) {
+                               fsgr->prepareVisualization();
+
+                               for(int k=0;k<=mLevel[simLev].lSizez-1;++k) {
+                                       for(int j=0;j<=mLevel[simLev].lSizey-1;++j) {
+                                               for(int i=isrc; i<mLevel[simLev].lSizex-1; i++) {
+                                                       //LbmFloat *cceldst = &QCELL(          msyLev ,idst+i,j,k, msySet,0);
+                                                       //LbmFloat *ccelsrc = &SIM_QCELL(fsgr, simLev ,isrc+i,j,k, simSet,0);
+                                                       //for(int l=0; l<dTotalNum; l++) {
+                                                               //RAC(cceldst,l) = RAC(ccelsrc,l);
+                                                       //}
+                                                       // both flag sets, make sure curr/other are same!?
+                                                       //RFLAG(msyLev ,idst+i,j,k, 0) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 0);
+                                                       //RFLAG(msyLev ,idst+i,j,k, 1) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 1);
+                                                       errMsg("TESTX"," "<<PRINT_VEC(idst+i,j,k)<<" < "<<PRINT_VEC(i,j,k) );
+                                                       *mpIso->lbmGetData(idst+i,j,k) = *simIso->lbmGetData(i,j,k);
+                                               }
+                                       }
+                               }
+
+                               } // simIso
+                       }
+               }
+       }
+#endif // LBM_INCLUDE_TESTSOLVERS==1
+
        return;
 }
 
@@ -315,7 +488,7 @@ int LbmFsgrSolver::initParticles() {
     //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid) ) 
     //&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid ) 
     //if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) { 
-    if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid) ) { 
+    if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid|CFUnused) ) { 
                        bool cellOk = true;
                        //? FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; }
                        if(!cellOk) continue;
@@ -426,10 +599,22 @@ int LbmFsgrSolver::initParticles() {
     /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \
                p->setType(newtype); 
 
+// tracer defines
+#define TRACE_JITTER 0.025
+#define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
+#define FFGET_NORM(var,dl) \
+                                                       if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFInter)){ (var) = QCELL_NB(lev,i,j,k,workSet,dl,dFfrac); } \
+                                                       else if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFFluid|CFUnused)){ (var) = 1.; } else (var) = 0.0;
+
+// float jitter
+#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
+#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) 
+
+
 void LbmFsgrSolver::advanceParticles() { 
   int workSet = mLevel[mMaxRefine].setCurr;
        LbmFloat vx=0.0,vy=0.0,vz=0.0;
-       LbmFloat rho, df[27]; //feq[27];
+       //LbmFloat df[27]; //feq[27];
 #define DEL_PART { \
        /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<"  ");  */ \
        p->setActive( false ); \
@@ -453,10 +638,6 @@ void LbmFsgrSolver::advanceParticles() {
        const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){ 
 
        // TODO use timestep size
-       //bool isIn,isOut,isInZ;
-       //const int cutval = 1+mCutoff/2; // TODO FIXME add half border!
-       //const int cutval = mCutoff/2; // TODO FIXME add half border!
-       //const int cutval = 0; // TODO FIXME add half border!
        const int cutval = mCutoff; // use full border!?
        int actCnt=0;
        if(this->mStepCnt%50==49) { mpParticles->cleanup(); }
@@ -518,27 +699,39 @@ void LbmFsgrSolver::advanceParticles() {
                    (p->getType()==PART_TRACER) ) {
 
                        // no interpol
-                       rho = vx = vy = vz = 0.0;
+                       vx = vy = vz = 0.0;
                        if(p->getStatus()&PART_IN) { // IN
                                if(k>=cutval) {
                                        if(k>this->mSizez-1-cutval) DEL_PART; 
-                                       FORDF0{
-                                               LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
-                                               df[l] = cdf;
-                                               rho += cdf; 
-                                               vx  += (this->dfDvecX[l]*cdf); 
-                                               vy  += (this->dfDvecY[l]*cdf);  
-                                               vz  += (this->dfDvecZ[l]*cdf);  
-                                       }
-
-                                       // remove gravity influence
-                                       const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
-                                       vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
-                                       vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
-                                       vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
 
-                                       if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid) ) {
+                                       if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFUnused) ) {
                                                // still ok
+                                               int partLev = mMaxRefine;
+                                               int si=i, sj=j, sk=k;
+                                               while(partLev>0 && RFLAG(partLev, si,sj,sk, workSet)&(CFUnused)) {
+                                                       partLev--;
+                                                       si/=2;
+                                                       sj/=2;
+                                                       sk/=2;
+                                               }
+                                               //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
+                                               LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet);
+                                               FORDF1{
+                                                       //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
+                                                       LbmFloat cdf = RAC(ccel, l);
+                                                       // TODO update below
+                                                       //df[l] = cdf;
+                                                       vx  += (this->dfDvecX[l]*cdf); 
+                                                       vy  += (this->dfDvecY[l]*cdf);  
+                                                       vz  += (this->dfDvecZ[l]*cdf);  
+                                               }
+
+                                               // remove gravity influence
+                                               const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
+                                               vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
+                                               vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
+                                               vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
+
                                        } else { // OUT
                                                // out of bounds, deactivate...
                                                // FIXME make fsgr treatment
@@ -548,12 +741,12 @@ void LbmFsgrSolver::advanceParticles() {
                                        // below 3d region, just rise
                                }
                        } else { // OUT
-#if LBM_INCLUDE_TESTSOLVERS==1
+#                              if LBM_INCLUDE_TESTSOLVERS==1
                                if(useff) { mpTest->handleParticle(p, i,j,k); }
                                else DEL_PART;
-#else // LBM_INCLUDE_TESTSOLVERS==1
+#                              else // LBM_INCLUDE_TESTSOLVERS==1
                                DEL_PART;
-#endif // LBM_INCLUDE_TESTSOLVERS==1
+#                              endif // LBM_INCLUDE_TESTSOLVERS==1
                                // TODO use x,y vel...?
                        }
 
@@ -577,15 +770,11 @@ void LbmFsgrSolver::advanceParticles() {
                                //LbmVec change = (fb+fd) *timestep / (pvolume*rhoAir)  *(timestep/cellsize);
                                //actCnt++; // should be after active test
                                if(actCnt<0) {
-                                       errMsg("\nPIT","BTEST1   vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
+                                       errMsg("PIT","BTEST1   vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
                                        errMsg("PIT","BTEST2        cellsize="<<cellsize<<" timestep="<<timestep<<" viscW="<<viscWater<<" ss/mb="<<(timestep/(pvolume*rhoAir)));
                                        errMsg("PIT","BTEST2        grav="<<rwgrav<<"  " );
                                        errMsg("PIT","BTEST2        change="<<(change)<<" fb="<<(fb)<<" fd="<<(fd)<<" ");
                                        errMsg("PIT","BTEST2        change="<<norm(change)<<" fb="<<norm(fb)<<" fd="<<norm(fd)<<" ");
-#if LOOPTEST==1
-                                       errMsg("PIT","BTEST2        n="<<n<<" "); // LOOPTEST! DEBUG
-#endif // LOOPTEST==1
-                                       errMsg("PIT","\n");
                                }
                                        
                                //v += change;
@@ -601,25 +790,15 @@ void LbmFsgrSolver::advanceParticles() {
                                if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true);
                                } else { p->setInFluid(false); }
 
-                               //if(( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) {
                                if( (( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) ||
                                                (( fflag&CFInter ) && (!(fflag&CFNoNbFluid)) ) ) {
                                        // only real fluid
-#define TRACE_JITTER 0.025
-#define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
-#define __TRACE_RAND 0.
-#if LBMDIM==3
+#                                      if LBMDIM==3
                                        p->advance( TRACE_RAND,TRACE_RAND,TRACE_RAND);
-#else
+#                                      else
                                        p->advance( TRACE_RAND,TRACE_RAND, 0.);
-#endif
+#                                      endif
 
-                               //} // test!?
-                               //if( fflag&(CFNoBndFluid) ) {
-                                       // ok
-                               //} else if( fflag&(CFEmpty|CFNoNbFluid) ) {
-                               //} else if( fflag&(CFEmpty) ) {
-                                       //v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
                                } else {
                                        // move inwards along normal, make sure normal is valid first
                                        const int lev = mMaxRefine;
@@ -629,25 +808,21 @@ void LbmFsgrSolver::advanceParticles() {
                                        if(i>=this->mSizex-1) { nx =  1.; nonorm = true; }
                                        if(j<=0)              { ny = -1.; nonorm = true; }
                                        if(j>=this->mSizey-1) { ny =  1.; nonorm = true; }
-#if LBMDIM==3
+#                                      if LBMDIM==3
                                        if(k<=0)              { nz = -1.; nonorm = true; }
                                        if(k>=this->mSizez-1) { nz =  1.; nonorm = true; }
-#endif // LBMDIM==3
-                       //mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
+#                                      endif // LBMDIM==3
                                        if(!nonorm) {
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dE,dFfrac); } else nv1 = 0.0;
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dW,dFfrac); } else nv2 = 0.0;
+                                               FFGET_NORM(nv1,dE); FFGET_NORM(nv2,dW);
                                                nx = 0.5* (nv2-nv1);
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dN,dFfrac); } else nv1 = 0.0;
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dS,dFfrac); } else nv2 = 0.0;
+                                               FFGET_NORM(nv1,dN); FFGET_NORM(nv2,dS);
                                                ny = 0.5* (nv2-nv1);
-#if LBMDIM==3
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dT,dFfrac); } else nv1 = 0.0;
-                                               if(RFLAG_NB(lev,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dB,dFfrac); } else nv2 = 0.0;
+#                                              if LBMDIM==3
+                                               FFGET_NORM(nv1,dT); FFGET_NORM(nv2,dB);
                                                nz = 0.5* (nv2-nv1);
-#else //LBMDIM==3
-                                               nz = 0.0;
-#endif //LBMDIM==3
+#                                              else // LBMDIM==3
+                                               nz = 0.;
+#                                              endif // LBMDIM==3
                                        } else {
                                                v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
                                        }
@@ -708,7 +883,7 @@ void LbmFsgrSolver::advanceParticles() {
                                        if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) {
                                                // still ok
                                        // shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
-                                       } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid) ){
+                                       } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFUnused) ){
                                                // FIXME make fsgr treatment
                                                P_CHANGETYPE(p, PART_FLOAT ); continue; 
                                                // jitter in cell to prevent stacking when hitting a steep surface
@@ -721,12 +896,12 @@ void LbmFsgrSolver::advanceParticles() {
                                        }
                                }
                        } else { // OUT
-#if LBM_INCLUDE_TESTSOLVERS==1
+#                              if LBM_INCLUDE_TESTSOLVERS==1
                                if(useff) { mpTest->handleParticle(p, i,j,k); }
                                else{ DEL_PART; }
-#else // LBM_INCLUDE_TESTSOLVERS==1
+#                              else // LBM_INCLUDE_TESTSOLVERS==1
                                { DEL_PART; }
-#endif // LBM_INCLUDE_TESTSOLVERS==1
+#                              endif // LBM_INCLUDE_TESTSOLVERS==1
                        }
 
                } // air particle
@@ -741,7 +916,7 @@ void LbmFsgrSolver::advanceParticles() {
 
                                if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
                                        // still ok
-                               } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
+                               } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
                        //errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
                                        //P_CHANGETYPE(p, PART_BUBBLE ); continue;
@@ -767,7 +942,7 @@ void LbmFsgrSolver::advanceParticles() {
                        //  test - delte on boundary!?
                        //if( (i<=cutval)||(i>=this->mSizex-1-cutval)|| (j<=cutval)||(j>=this->mSizey-1-cutval)) { DEL_PART; } // DEBUG TEST
 
-#if LBM_INCLUDE_TESTSOLVERS==1
+#                              if LBM_INCLUDE_TESTSOLVERS==1
                        /*LbmFloat prob = 1.0;
                        // vanishing
                        prob = (rand()/(RAND_MAX+1.0));
@@ -783,30 +958,26 @@ void LbmFsgrSolver::advanceParticles() {
                        //? if(k>this->mSizez*3/4) {    if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
                        // */
 
-#else // LBM_INCLUDE_TESTSOLVERS==1
-#endif // LBM_INCLUDE_TESTSOLVERS==1
+#                              endif // LBM_INCLUDE_TESTSOLVERS==1
 
                        if(p->getStatus()&PART_IN) { // IN
                                //if((k<cutval)||(k>this->mSizez-1-cutval)) DEL_PART; 
                                if(k<cutval) DEL_PART; 
                                if(k>this->mSizez-1-cutval) { P_CHANGETYPE(p, PART_DROP ); continue; } // create new drop
                                //ntlVec3Gfx v = getVelocityAt(i,j,k);
-                               rho = vx = vy = vz = 0.0;
+                               vx = vy = vz = 0.0;
 
                                // define from particletracer.h
 #if MOVE_FLOATS==1
                                const int DEPTH_AVG=3; // only average interface vels
                                int ccnt=0;
                                for(int kk=0;kk<DEPTH_AVG;kk+=1) {
-                               //for(int kk=1;kk<DEPTH_AVG;kk+=1) {
                                        if((k-kk)<1) continue;
-                                       //if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue;
                                        if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue;
                                        ccnt++;
-                                       FORDF0{
+                                       FORDF1{
                                                LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l);
-                                               df[l] = cdf;
-                                               //rho += cdf; 
+                                               //df[l] = cdf;
                                                vx  += (this->dfDvecX[l]*cdf); 
                                                vy  += (this->dfDvecY[l]*cdf);  
                                                vz  += (this->dfDvecZ[l]*cdf);  
@@ -824,7 +995,7 @@ void LbmFsgrSolver::advanceParticles() {
                                vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
 
                                //bool delfloat = false;
-                               if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
+                               if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
                                        // in fluid cell
                                        /*if((1) && (k<this->mSizez-3) && 
                                                        (
@@ -888,16 +1059,10 @@ void LbmFsgrSolver::advanceParticles() {
                                // use half butoff border 1/8
                                int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5);
                                if(maxdw<3) maxdw=3;
-#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
-#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) 
-                               //if(ABS(i-(               cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
-                               //if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
                                if((j>=0)&&(j<=this->mSizey-1)) {
                                        if(ABS(i-(               cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(               cutval))), 0.,0.); }
                                        if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); }
                                }
-//#undef FLOAT_JITTER_BND
-#undef FLOAT_JITTBNDRAND
                                //if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval)
                        }
                }  // PART_FLOAT
@@ -1306,11 +1471,15 @@ int LbmFsgrSolver::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((ii==-1)&&(ij==0)) {
+               // special case for main loop, ok
+       } else {
+               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(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<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; }
@@ -1338,11 +1507,15 @@ int LbmFsgrSolver::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((ii==-1)&&(ij==0)) {
+               // special case for main loop, ok
+       } else {
+               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(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<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; }
@@ -1561,7 +1734,10 @@ void LbmFsgrSolver::debugDisplayNode(int dispset, CellIdentifierInterface* cell
 void LbmFsgrSolver::lbmDebugDisplay(int dispset) {
        // DEBUG always display testdata
 #if LBM_INCLUDE_TESTSOLVERS==1
-       if(mUseTestdata){ mpTest->testDebugDisplay(dispset); }
+       if(mUseTestdata){ 
+               cpDebugDisplay(dispset); 
+               mpTest->testDebugDisplay(dispset); 
+       }
 #endif // LBM_INCLUDE_TESTSOLVERS==1
        if(dispset<=FLUIDDISPNothing) return;
        //if(!dispset->on) return;