solver_main.cpp

Go to the documentation of this file.

/******************************************************************************
 *
 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
 * Copyright 2003-2006 Nils Thuerey
 *
 * Standard LBM Factory implementation
 *
 *****************************************************************************/

#include "solver_class.h"
#include "solver_relax.h"
#include "particletracer.h"
#include "loop_tools.h"
#include <stdlib.h>

/*****************************************************************************/
/*****************************************************************************/

double globdfcnt;
double globdfavg[19];
double globdfmax[19];
double globdfmin[19];
extern int glob_mpindex,glob_mpnum;
extern bool globOutstrForce;

// simulation object interface
void LbmFsgrSolver::step() { 
    stepMain();
}

// lbm main step
void messageOutputForce(string from);
void LbmFsgrSolver::stepMain() { 
    myTime_t timestart = getTime();

    initLevelOmegas();
    markedClearList(); // DMC clearMarkedCellsList

    // safety check, counter reset
    mNumUsedCells = 0;
    mNumInterdCells = 0;
    mNumInvIfCells = 0;

  //debugOutNnl("LbmFsgrSolver::step : "<<mStepCnt, 10);
  if(!mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime, 10); }
    //debMsgDirect(  "LbmFsgrSolver::step : "<<mStepCnt<<" ");
    //myTime_t timestart = 0;
    //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 

    // time adapt
    mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;

    // init moving bc's, can change mMaxVlen
    initMovingObstacles(false);
    
    // handle fluid control 
    handleCpdata();

    // important - keep for tadap
    LbmFloat lastMass = mCurrentMass;
    mCurrentMass = mFixMass; // reset here for next step
    mCurrentVolume = 0.0;
    
    //change to single step advance!
    int levsteps = 0;
    int dsbits = mStepCnt ^ (mStepCnt-1);
    //errMsg("S"," step:"<<mStepCnt<<" s-1:"<<(mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
    for(int lev=0; lev<=mMaxRefine; lev++) {
        //if(! (mStepCnt&(1<<lev)) ) {
        if( dsbits & (1<<(mMaxRefine-lev)) ) {
            //errMsg("S"," l"<<lev);

            if(lev==mMaxRefine) {
                // always advance fine level...
                fineAdvance();
            } else {
                adaptGrid(lev);
                coarseRestrictFromFine(lev);
                coarseAdvance(lev);
            }
#if FSGR_OMEGA_DEBUG==1
            errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
            mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
#endif // FSGR_OMEGA_DEBUG==1
            levsteps++;
        }
        mCurrentMass   += mLevel[lev].lmass;
        mCurrentVolume += mLevel[lev].lvolume;
    }

  // prepare next step
    mStepCnt++;


    // some dbugging output follows
    // calculate MLSUPS
    myTime_t timeend = getTime();

    mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
    mAvgNumUsedCells += mNumUsedCells;
    mMLSUPS = ((double)mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000.0);
    if(mMLSUPS>10000){ mMLSUPS = -1; }
    //else { mAvgMLSUPS += mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
    
    LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
    if(totMLSUPS>10000) totMLSUPS = -1;
    mNumInvIfTotal += mNumInvIfCells; // debug

  // do some formatting 
  if(!mSilent){ 
        int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt);
    debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime<<
            " cur-mlsups:"<<mMLSUPS<< //" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< 
            " totcls:"<<mNumUsedCells<< " avgcls:"<< avgcls<< 
            " intd:"<<mNumInterdCells<< " invif:"<<mNumInvIfCells<< 
            " invift:"<<mNumInvIfTotal<< " fsgrcs:"<<mNumFsgrChanges<< 
            " filled:"<<mNumFilledCells<<", emptied:"<<mNumEmptiedCells<< 
            " mMxv:"<<PRINT_VEC(mMxvx,mMxvy,mMxvz)<<", tscnts:"<<mTimeSwitchCounts<< 
            //" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
            " probs:"<<mNumProblems<< " simt:"<<mSimulationTime<< 
            " took:"<< getTimeString(timeend-timestart)<<
            " for '"<<mName<<"' " , 10);
    } else { debMsgDirect("."); }

    if(mStepCnt==1) {
        mMinNoCells = mMaxNoCells = mNumUsedCells;
    } else {
        if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells;
        if(mNumUsedCells<mMinNoCells) mMinNoCells = mNumUsedCells;
    }
    
    // mass scale test
    if((mMaxRefine>0)&&(mInitialMass>0.0)) {
        LbmFloat mscale = mInitialMass/mCurrentMass;

        mscale = 1.0;
        const LbmFloat dchh = 0.001;
        if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
        if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;

        // use mass rescaling?
        // with float precision this seems to be nonsense...
        const bool MREnable = false;

        const int MSInter = 2;
        static int mscount = 0;
        if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
            // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
            // mass rescale MASS RESCALE check
            errMsg("MDTDD","\n\n");
            errMsg("MDTDD","FORCE RESCALE MASS! "
                    <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
                    <<" step:"<<mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
                    );
            errMsg("MDTDD","\n\n");

            mscount++;
            for(int lev=mMaxRefine; lev>=0 ; lev--) {
                //for(int workSet = 0; workSet<=1; workSet++) {
                int wss = 0;
                int wse = 1;
#if COMPRESSGRIDS==1
                if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
#endif // COMPRESSGRIDS==1
                for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT

                    FSGR_FORIJK1(lev) {
                        if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
                            ) {

                            FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
                            QCELL(lev, i,j,k,workSet, dMass) *= mscale;
                            QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;

                        } else {
                            continue;
                        }
                    }
                }
                mLevel[lev].lmass *= mscale;
            }
        } 

        mCurrentMass *= mscale;
    }// if mass scale test */
    else {
        // use current mass after full step for initial setting
        if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
            mInitialMass = mCurrentMass;
            debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
        }
    }

#if LBM_INCLUDE_TESTSOLVERS==1
    if((mUseTestdata)&&(mInitDone)) { handleTestdata(); }
    mrExchange();
#endif

    // advance positions with current grid
    advanceParticles();
    if(mpParticles) {
        mpParticles->checkDumpTextPositions(mSimulationTime);
        mpParticles->checkTrails(mSimulationTime);
    }

    // one of the last things to do - adapt timestep
    // was in fineAdvance before... 
    if(mTimeAdap) {
        adaptTimestep();
    } // time adaptivity


#ifndef WIN32
    // good indicator for instabilities
    if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
    if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
#endif // WIN32

    // output total step time
    myTime_t timeend2 = getTime();
    double mdelta = (lastMass-mCurrentMass);
    if(ABS(mdelta)<1e-12) mdelta=0.;
    double effMLSUPS = ((double)mNumUsedCells / ((timeend2-timestart)/(double)1000.0) ) / (1000000.0);
    if(mInitDone) {
        if(effMLSUPS>10000){ effMLSUPS = -1; }
        else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
    }
    
    debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<<glob_mpindex<<" step:"<<mStepCnt
            <<" dccd="<< mCurrentMass
            //<<",d"<<mdelta
            //<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
            //<<"/"<<mCurrentVolume<<"(fix="<<mFixMass<<",ini="<<mInitialMass<<"), "
            <<" effMLSUPS=("<< effMLSUPS
            <<",avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "
            <<" took totst:"<< getTimeString(timeend2-timestart), 3);
    // nicer output
    //debMsgDirect(std::endl); 
    // */

    messageOutputForce("");
 //#endif // ELBEEM_PLUGIN!=1
}

#define NEWDEBCHECK(str) \
    if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \
        if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \
        for(int l=0;l<dTotalNum;l++) { \
            if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\
        }/*for*/ \
        }/*if*/ \
    } }

void LbmFsgrSolver::fineAdvance()
{
    // do the real thing...
    //NEWDEBCHECK("t1");
    mainLoop( mMaxRefine );
    if(mUpdateFVHeight) {
        // warning assume -Y gravity...
        mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
        if(mFVHeight<1.0) mFVHeight = 1.0;
        mpParam->setFluidVolumeHeight(mFVHeight);
    }

    // advance time before timestep change
    mSimulationTime += mpParam->getTimestep();
    // time adaptivity
    mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
    //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
    if(!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;
  mLevel[mMaxRefine].setCurr   ^= 1;
  mLevel[mMaxRefine].lsteps++;

    // flag init... (work on current set, to simplify flag checks)
    reinitFlags( mLevel[mMaxRefine].setCurr );
    if(!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 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 (mSizez/4)
#else // LBMDIM==2
// off
#define SLOWDOWNREGION 10 
#endif // LBMDIM==2
#define P_LCSMQO 0.01

/*****************************************************************************/
/*****************************************************************************/
void 
LbmFsgrSolver::mainLoop(int lev)
{
    // loops over _only inner_ cells  -----------------------------------------------------------------------------------
    
    // 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->mFarfMode>0)) { return; }
#   endif // ELBEEM_PLUGIN!=1
    
  // main loop region
    const bool doReduce = true;
    const int gridLoopBound=1;
    GRID_REGION_INIT();
#if PARALLEL==1
#pragma omp parallel default(shared) \
  reduction(+: \
      calcCurrentMass,calcCurrentVolume, \
        calcCellsFilled,calcCellsEmptied, \
        calcNumUsedCells )
    GRID_REGION_START();
#else // PARALLEL==1
    GRID_REGION_START();
#endif // PARALLEL==1

    // local to main
    CellFlagType nbflag[LBM_DFNUM]; 
    int oldFlag, newFlag, nbored;
    LbmFloat m[LBM_DFNUM];
    LbmFloat rho, ux, uy, uz, tmp, usqr;

    // smago vars
    LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
    
    // ifempty cell conversion flags
    bool iffilled, ifemptied;
    LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
    int recons[LBM_DFNUM];   // reconstruct this DF?
    int numRecons;           // how many are reconstructed?

    LbmFloat mass=0., change=0., lcsmqo=0.;
    rho= ux= uy= uz= usqr= tmp= 0.; 
    lcsmqadd = lcsmomega = 0.;
    FORDF0{ lcsmeq[l] = 0.; }

    // ---
    // now stream etc.
    // use //template functions for 2D/3D

    GRID_LOOP_START();
        // loop start
        // stream from current set to other, then collide and store
        //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);

#       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) ) {
            errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
                    RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
                    (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
          (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
          (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
                    ); 
            CAUSE_PANIC;
        }   
        if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
            (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
            errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
          (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
          (long)(ccel)<<","<<(long)(tcel)<<" "
                    ); 
            CAUSE_PANIC;
        }   
#       endif
        oldFlag = *pFlagSrc;
        
        // old INTCFCOARSETEST==1
        if( (oldFlag & (CFGrFromCoarse)) ) { 
            if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
                FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
            } else {
                interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
                calcNumUsedCells++;
            }
            continue; // interpolateFineFromCoarse test!
        } // interpolateFineFromCoarse test! 
    
        if(oldFlag & (CFMbndInflow)) {
            // fluid & if are ok, fill if later on
            int isValid = oldFlag & (CFFluid|CFInter);
            const LbmFloat iniRho = 1.0;
            const int OId = oldFlag>>24;
            if(!isValid) {
                // make new if cell
                const LbmVec vel(mObjectSpeeds[OId]);
                // TODO add OPT3D treatment
                FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
                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++;
                mInitialMass += iniRho;
                // dont treat cell until next step
                continue;
            } 
        } 
        else  // these are exclusive
        if(oldFlag & (CFMbndOutflow)) {
            int isnotValid = oldFlag & (CFFluid);
            if(isnotValid) {
                // remove fluid cells, shouldnt be here anyway
                LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
                mInitialMass -= fluidRho;
                const LbmFloat iniRho = 0.0;
                RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
                RAC(tcel, dFlux) = FLUX_INIT;
                changeFlag(lev, i,j,k, TSET(lev), CFInter);

                // same as ifemptied for if below
                LbmPoint oemptyp; oemptyp.flag = 0;
                oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
                LIST_EMPTY(oemptyp);
                calcCellsEmptied++;
                continue;
            }
        }

        if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
            *pFlagDst = oldFlag;
            continue;
        }
        /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
            OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
            RAC(tcel,dFfrac) = 1.0; 
            *pFlagDst = (CellFlagType)oldFlag; // newFlag;
            calcCurrentMass += rho; calcCurrentVolume += 1.0;
            calcNumUsedCells++;
            continue;
        }// TEST ME FASTER? */

        // only neighbor flags! not own flag
        nbored = 0;
        
#if OPT3D==0
        FORDF1 {
            nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
            nbored |= nbflag[l];
        } 
#else
        nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
        nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
        nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
        nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
        nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];

        nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
        nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
        nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];

        nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
        nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];

        nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
      nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
        nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];

        nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
        nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
        nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
        nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
        nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
        // */
#endif

        // pointer to destination cell
        calcNumUsedCells++;

        // FLUID cells 
        if( oldFlag & CFFluid ) { 
            // only standard fluid cells (with nothing except fluid as nbs

            if(oldFlag&CFMbndInflow) {
                // force velocity for inflow, necessary to have constant direction of flow
                // FIXME , test also set interface cells?
                const int OId = oldFlag>>24;
                //? DEFAULT_STREAM;
                //const LbmFloat fluidRho = 1.0;
                // for submerged inflows, streaming would have to be performed...
                LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
                const LbmVec vel(mObjectSpeeds[OId]);
                ux=vel[0], uy=vel[1], uz=vel[2]; 
                usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
                FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
                rho = fluidRho;
                //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
            } else {
                if(nbored&CFBnd) {
                    DEFAULT_STREAM;
                    //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
                    DEFAULT_COLLIDEG(mLevel[lev].gravity);
                    oldFlag &= (~CFNoBndFluid);
                } else {
                    // do standard stream/collide
                    OPTIMIZED_STREAMCOLLIDE;
                    oldFlag |= CFNoBndFluid;
                } 
            }

            PERFORM_USQRMAXCHECK;
            // "normal" fluid cells
            RAC(tcel,dFfrac) = 1.0; 
            *pFlagDst = (CellFlagType)oldFlag; // newFlag;
            calcCurrentMass += rho; 
            calcCurrentVolume += 1.0;
            continue;
        }
        
        newFlag  = oldFlag;
        // make sure here: always check which flags to really unset...!
        newFlag = newFlag & (~( 
                    CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
                    | CFNoInterpolSrc
                    | CFNoBndFluid
                    ));
        if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
            newFlag |= CFNoBndFluid;
        }
        /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
            // explicitly check for noslip neighbors
            bool hasnoslipnb = false;
            FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
            if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
        } // */

        // store own dfs and mass
        mass = RAC(ccel,dMass);

        // WARNING - only interface cells arrive here!
        // read distribution funtions of adjacent cells = stream step
        DEFAULT_STREAM;

        if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
        if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }

        // calculate mass exchange for interface cells 
        LbmFloat myfrac = RAC(ccel,dFfrac);
        if(myfrac<0.) myfrac=0.; // NEWSURFT
#       define nbdf(l) m[ this->dfInv[(l)] ]

        // update mass 
        // only do boundaries for fluid cells, and interface cells without
        // any fluid neighbors (assume that interface cells _with_ fluid
        // neighbors are affected enough by these) 
        // which Df's have to be reconstructed? 
        // for fluid cells - just the f_i difference from streaming to empty cells  ----
        numRecons = 0;
        bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
        //onlyBndnb = false; // DEBUG test off

        FORDF1 { // dfl loop
            recons[l] = 0;
            nbfracs[l] = 0.0;
            // finally, "normal" interface cells ----
            if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
                change = nbdf(l) - MYDF(l);
            }
            // interface cells - distuingish cells that shouldn't fill/empty 
            else if( nbflag[l] & CFInter ) {
                
                LbmFloat mynbfac,nbnbfac;
                // NEW TEST t1
                // t2 -> off
                if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
                    mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
                    nbnbfac = 1.0/mynbfac;
                    onlyBndnb = false;
                } else {
                    mynbfac = nbnbfac = 1.0; // switch calc flux off
                    goto changeDefault;  // NEWSURFT
                    //change = 0.; goto changeDone;  // NEWSURFT
                }
                //mynbfac = nbnbfac = 1.0; // switch calc flux off t3

                // perform interface case handling
                if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
                switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
                    case 0: 
                        // we are a normal cell so... 
                        switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
                            case CFNoNbFluid: 
                                // just fill current cell = empty neighbor 
                                change = nbnbfac*nbdf(l) ; goto changeDone; 
                            case CFNoNbEmpty: 
                                // just empty current cell = fill neighbor 
                                change = - mynbfac*MYDF(l) ; goto changeDone; 
                        }
                        break;

                    case CFNoNbFluid: 
                        // we dont have fluid nb's so...
                        switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
                            case 0: 
                            case CFNoNbEmpty: 
                                // we have no fluid nb's -> just empty
                                change = - mynbfac*MYDF(l) ; goto changeDone; 
                        }
                        break;

                    case CFNoNbEmpty: 
                        // we dont have empty nb's so...
                        switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
                            case 0: 
                            case CFNoNbFluid: 
                                // we have no empty nb's -> just fill
                                change = nbnbfac*nbdf(l); goto changeDone; 
                        }
                        break;
                }} // inter-inter exchange

            changeDefault: ;
                // just do normal mass exchange...
                change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
            changeDone: ;
                nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
                if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
                change *=  (myfrac + nbfracs[l]) * 0.5;
            } // the other cell is interface

            // last alternative - reconstruction in this direction
            else {
                // empty + bnd case
                recons[l] = 1; 
                numRecons++;
                change = 0.0; 
            }

            // modify mass at SRCS
            mass += change;
        } // l
        // normal interface, no if empty/fluid

        // computenormal
        LbmFloat surfaceNormal[3];
        computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal);

        if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) {
            // normal ok and usable...
            FORDF1 {
                if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2])  // dot Dvec,norml
                        > LBM_EPSILON) {
                    recons[l] = 2; 
                    numRecons++;
                } 
            }
        }

        // calculate macroscopic cell values
        LbmFloat oldUx, oldUy, oldUz;
        LbmFloat oldRho; // OLD rho = ccel->rho;
#       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 );
            }
        }
        ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
        usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // 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 )
                + RAC(ccel,dB ) + RAC(ccel,dNE)
                + RAC(ccel,dNW) + RAC(ccel,dSE)
                + RAC(ccel,dSW) + RAC(ccel,dNT)
                + RAC(ccel,dNB) + RAC(ccel,dST)
                + RAC(ccel,dSB) + RAC(ccel,dET)
                + RAC(ccel,dEB) + RAC(ccel,dWT)
                + RAC(ccel,dWB);

        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)
                + 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)
                + RAC(ccel,dNT) - RAC(ccel,dNB)
                + RAC(ccel,dST) - RAC(ccel,dSB)
                + RAC(ccel,dET) - RAC(ccel,dEB)
                + RAC(ccel,dWT) - RAC(ccel,dWB);

        // now reconstruction
        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
        if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
        if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
        if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
        if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
        if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
        if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
        if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
        if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
        if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
        if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
        if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
        if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
        if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
        if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
        if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
        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       


        // inflow bc handling
        if(oldFlag & (CFMbndInflow)) {
            // fill if cells in inflow region
            if(myfrac<0.5) { 
                mass += 0.25; 
                mInitialMass += 0.25;
            }
            const int OId = oldFlag>>24;
            const LbmVec vel(mObjectSpeeds[OId]);
            ux=vel[0], uy=vel[1], uz=vel[2]; 
            //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
            //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
            rho = REFERENCE_PRESSURE;
            FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
            //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
        }  else { 
            // NEWSURFT, todo optimize!
            if(onlyBndnb) { //if(usqr<0.001*0.001) {
                rho = ux = uy = uz = 0.;
                FORDF0{ 
                    rho += m[l]; 
                    ux  += (this->dfDvecX[l]*m[l]); 
                    uy  += (this->dfDvecY[l]*m[l]);  
                    uz  += (this->dfDvecZ[l]*m[l]);  
                }
                FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
            } else {// NEWSURFT */
                if(usqr>0.3*0.3) { 
                    // force reset! , warning causes distortions...
                    FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); }
                } else {
                // normal collide
                // mass streaming done... do normal collide
                LbmVec grav = mLevel[lev].gravity*mass;
                DEFAULT_COLLIDEG(grav);
                PERFORM_USQRMAXCHECK; }
                // rho init from default collide necessary for fill/empty check below
            } // test
        }

        // testing..., particle generation
        // also check oldFlag for CFNoNbFluid, test
        // for inflow no pargen test // NOBUBBB!
        if((mInitDone) 
                // dont allow new if cells, or submerged ones
                && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) 
                // dont try to subtract from empty cells
                && (mass>0.) && (mPartGenProb>0.0)) {
            bool doAdd = true;
            bool bndOk=true;
            if( (i<cutMin)||(i>mSizex-cutMin)||
                    (j<cutMin)||(j>mSizey-cutMin)||
                    (k<cutMin)||(k>mSizez-cutMin) ) { bndOk=false; }
            if(!bndOk) doAdd=false;
            
            LbmFloat prob = (rand()/(RAND_MAX+1.0));
            LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333;

            // physical drop model
            if(mPartUsePhysModel) {
                LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
                LbmFloat rux = (ux * realWorldFac);
                LbmFloat ruy = (uy * realWorldFac);
                LbmFloat ruz = (uz * realWorldFac);
                LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
                basethresh *= rl;

                // 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.;

                //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
                if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
                    // add
                } else {
                    doAdd = false; // dont...
                }

                // "wind" disturbance
                // use realworld relative velocity here instead?
                if( (doAdd && 
                        ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
                        ||(k>mSizez-SLOWDOWNREGION)   ) {
                    LbmFloat nuz = uz;
                    if(k>mSizez-SLOWDOWNREGION) {
                        // special case
                        LbmFloat zfac = (LbmFloat)( k-(mSizez-SLOWDOWNREGION) );
                        zfac /= (LbmFloat)(SLOWDOWNREGION);
                        nuz += (1.0) * zfac; // check max speed? OFF?
                        //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac);
                    } else {
                        // normal probability
                        //? LbmFloat fac = P_LCSMQO-lcsmqo;
                        //? jdf *= fac;
                    }
                    FORDF1 {
                        const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0));
                        // TODO  use wind velocity?
                        if(jdf>0.025) {
                        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 "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
                } // wind disturbance
            } // mPartUsePhysModel
            else {
                // empirical model
                //if((prob<basethresh) && (lcsmqo>0.0095)) { // add
                if((prob<basethresh) && (lcsmqo>0.012)) { // add
                } else { doAdd = false; }// dont...
            } 


            // remove noise
            if(usqr<0.0001) doAdd=false;   // TODO check!?

            // dont try to subtract from empty cells
            // ensure cell has enough mass for new drop
            LbmFloat newPartsize = 1.0;
            if(mPartUsePhysModel) {
                // 1-10
                newPartsize += 9.0* (rand()/(RAND_MAX+1.0));
            } else {
                // 1-5, overall size has to be less than
                // .62 (ca. 0.5) to make drops significantly smaller 
                // than a full cell!
                newPartsize += 4.0* (rand()/(RAND_MAX+1.0));
            }
            LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); //PARTMASS(mPartDropMassSub*newPartsize); // mass: 4/3 pi r^3 rho
            while(dropmass>mass) {
                newPartsize -= 0.2;
                dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize);
            }
            if(newPartsize<=1.) doAdd=false;

            if( (doAdd)  ) { // init new particle
                for(int s=0; s<1; s++) { // one part!
                const LbmFloat posjitter = 0.05;
                const LbmFloat posjitteroffs = posjitter*-0.5;
                LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
                LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
                LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));

                const LbmFloat jitterstr = 1.0;
                const LbmFloat jitteroffs = jitterstr*-0.5;
                LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
                LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
                LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));

                // average normal & velocity 
                // -> mostly along velocity dir, many into surface
                // fluid velocity (not normalized!)
                LbmVec flvelVel = LbmVec(ux,uy,uz);
                LbmFloat flvelLen = norm(flvelVel);
                // surface normal
                LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]);
                normalize(normVel);
                LbmFloat normScale = (0.01+flvelLen);
                // jitter vector, 0.2 * flvel
                LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1;
                // weighten velocities
                const LbmFloat flvelWeight = 0.9;
                LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel; 

                // offset towards surface (hide popping)
                jpx += -normVel[0]*0.4;
                jpy += -normVel[1]*0.4;
                jpz += -normVel[2]*0.4;

                LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz;
                int type=0;
                type = PART_DROP;

#               if LBMDIM==2
                newpartVel[2]=0.; srck=0.5;
#               endif // LBMDIM==2
                // subtract drop mass
                mass -= dropmass;
                // init new particle
                {
                    ParticleObject np( ntlVec3Gfx(srci,srcj,srck) );
                    np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]);
                    np.setStatus(PART_IN);
                    np.setType(type);
                    //if((s%3)==2) np.setType(PART_FLOAT);
                    np.setSize(newPartsize);
                    //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize );
                    //errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel );
                    FSGR_ADDPART(np);
                } // new part
            //errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel()  );
                } // multiple parts
            } // doAdd
        } // */


        // interface cell filled or emptied?
        iffilled = ifemptied = 0;
        // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
        // interface cell if full?
        if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
        // interface cell if empty?
        if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }

        if(oldFlag & (CFMbndOutflow)) {
            mInitialMass -= mass;
            mass = myfrac = 0.0;
            iffilled = 0; ifemptied = 1;
        }

        // looks much nicer... LISTTRICK
#       if FSGR_LISTTRICK==1
        //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
        if(newFlag&CFNoBndFluid) { // test NEW TEST
            if(!iffilled) {
                // remove cells independent from amount of change...
                if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
                        ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
                    //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
                    iffilled = 1; 
                } 
            }
            if(!ifemptied) {
                if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
                        ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
                    //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
                    ifemptied = 1; 
                } 
            }
        } // nobndfluid only */
#       endif
        //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
        

        // now that all dfs are known, handle last changes
        if(iffilled) {
            LbmPoint filledp; filledp.flag=0;
            if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
            filledp.x = i; filledp.y = j; filledp.z = k;
            LIST_FULL(filledp);
            //mNumFilledCells++; // DEBUG
            calcCellsFilled++;
        }
        else if(ifemptied) {
            LbmPoint emptyp; emptyp.flag=0;
            if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
            emptyp.x = i; emptyp.y = j; emptyp.z = k;
            LIST_EMPTY(emptyp);
            //mNumEmptiedCells++; // DEBUG
            calcCellsEmptied++;
        } 
        // dont cutoff values -> better cell conversions
        RAC(tcel,dFfrac)   = (mass/rho);

        // init new flux value
        float flux = FLUX_INIT; // dxqn on
        if(newFlag&CFNoBndFluid) {
            //flux = 50.0; // extreme on
            for(int nn=1; nn<this->cDfNum; nn++) { 
                if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
            }
            // optical hack - smooth slow moving
            // surface regions
            if(usqr< sssUsqrLimit) {
            for(int nn=1; nn<this->cDfNum; nn++) { 
                if(nbfracs[nn]!=0.0) {
                    LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
                    if(curSmooth>1.0) curSmooth=1.0;
                    flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
                }
            } }
            // NEW TEST */
        }
        // flux = FLUX_INIT; // calc flux off
        QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */

        // perform mass exchange with streamed values
        QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
        // set new flag 
        *pFlagDst = (CellFlagType)newFlag;
        calcCurrentMass += mass; 
        calcCurrentVolume += RAC(tcel,dFfrac);

        // interface cell handling done...

#if PARALLEL!=1
    GRID_LOOPREG_END();
#else // PARALLEL==1
#include "paraloopend.h" // = GRID_LOOPREG_END();
#endif // PARALLEL==1

    // write vars from computations to class
    mLevel[lev].lmass    = calcCurrentMass;
    mLevel[lev].lvolume  = calcCurrentVolume;
    mNumFilledCells  = calcCellsFilled;
    mNumEmptiedCells = calcCellsEmptied;
    mNumUsedCells = calcNumUsedCells;
}



void 
LbmFsgrSolver::preinitGrids()
{
    const int lev = mMaxRefine;
    const bool doReduce = false;
    const int gridLoopBound=0;

    // preinit both grids
    for(int s=0; s<2; s++) {
    
        GRID_REGION_INIT();
#if PARALLEL==1
#pragma omp parallel default(shared) \
  reduction(+: \
      calcCurrentMass,calcCurrentVolume, \
        calcCellsFilled,calcCellsEmptied, \
        calcNumUsedCells )
#endif // PARALLEL==1
        GRID_REGION_START();
        GRID_LOOP_START();
            for(int l=0; l<dTotalNum; l++) { 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

        /* 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
#pragma omp parallel default(shared) \
  reduction(+: \
      calcCurrentMass,calcCurrentVolume, \
        calcCellsFilled,calcCellsEmptied, \
        calcNumUsedCells )
#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();
        //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
        const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
        if( (currFlag & (CFEmpty|CFBnd)) ) continue;

        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;
            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

    /* dummy remove warnings */ 
    calcCurrentMass = calcCurrentVolume = 0.;
    calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
    
    // change grid
  mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
  mLevel[mMaxRefine].setCurr   ^= 1;
}


/******************************************************************************
 * work on lists from updateCellMass to reinit cell flags
 *****************************************************************************/

LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
    //return 0.0; // test
    int level = mMaxRefine;
    LbmFloat     *ccel  = RACPNT(level, i,j,k, workSet);

    // computenormal
    CellFlagType *cflag = &RFLAG(level, i,j,k, workSet);
    LbmFloat n[3];
    computeFluidSurfaceNormal(ccel,cflag, n);
    LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2];

    LbmFloat ret = 1.0;
    // forward direction, add mass (for filling cells):
    if(dirForw) {
        if(scal<LBM_EPSILON) ret = 0.0;
        else ret = scal;
    } else {
        // backward for emptying
        if(scal>-LBM_EPSILON) ret = 0.0;
        else ret = scal * -1.0;
    }
    //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
    return ret;
}

// warning - normal compuations are without
//   boundary checks &
//   normalization
void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) {
    const int level = mMaxRefine;
    LbmFloat nx,ny,nz, nv1,nv2;
    const CellFlagType flagE = *(cflagpnt+1);
    const CellFlagType flagW = *(cflagpnt-1);
    if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } 
    else if(flagE &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } 
    else if(flagW &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    nx = 0.5* (nv2-nv1);

    const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
    const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
    if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
    else if(flagN &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
    else if(flagS &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    ny = 0.5* (nv2-nv1);

#if LBMDIM==3
    const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
    const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
    if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
    else if(flagT &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
    else if(flagB &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    nz = 0.5* (nv2-nv1);
#else //LBMDIM==3
    nz = 0.0;
#endif //LBMDIM==3

    // return vals
    snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
    LbmFloat nx=0.,ny=0.,nz=0.;
    ccel = NULL; cflagpnt=NULL; // remove warning
    snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
    const int level = mMaxRefine;
    LbmFloat nx,ny,nz, nv1,nv2;
    ccel = NULL; // remove warning

    const CellFlagType flagE = *(cflagpnt+1);
    const CellFlagType flagW = *(cflagpnt-1);
    if(flagE &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagW &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    nx = 0.5* (nv2-nv1);

    const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
    const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
    if(flagN &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagS &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    ny = 0.5* (nv2-nv1);

#if LBMDIM==3
    const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
    const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
    if(flagT &(CFBnd)){ nv1 = 1.; }
    else nv1 = 0.0;
    if(flagB &(CFBnd)){ nv2 = 1.; }
    else nv2 = 0.0;
    nz = 0.5* (nv2-nv1);
#else //LBMDIM==3
    nz = 0.0;
#endif //LBMDIM==3

    // return vals
    snret[0]=nx; snret[1]=ny; snret[2]=nz;
}
void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) {
    bool nonorm = false;
    LbmFloat nx=0.,ny=0.,nz=0.;
    if(i<=0)        { nx =  1.; nonorm = true; }
    if(i>=mSizex-1) { nx = -1.; nonorm = true; }
    if(j<=0)        { ny =  1.; nonorm = true; }
    if(j>=mSizey-1) { ny = -1.; nonorm = true; }
#   if LBMDIM==3
    if(k<=0)        { nz =  1.; nonorm = true; }
    if(k>=mSizez-1) { nz = -1.; nonorm = true; }
#   endif // LBMDIM==3
    if(!nonorm) {
        // in domain, revert to helper, use setCurr&mMaxRefine
        LbmVec bnormal;
        CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
        LbmFloat     *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
        computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]);
        // TODO check if there is a normal near here?
        // use wider range otherwise...
        snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0];
        return;
    }
    snret[0]=nx; snret[1]=ny; snret[2]=nz;
}

void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
#if FSGR_STRICT_DEBUG==10
    // dangerous, this can change the simulation...
  /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
       iter != mListNewInter.end(); iter++ ) {
    if(ni!=iter->x) continue;
    if(nj!=iter->y) continue;
    if(nk!=iter->z) continue;
        // all 3 values match... skip point
        return;
    } */
#endif // FSGR_STRICT_DEBUG==1
    // store point
    LbmPoint newinter; newinter.flag = 0;
    newinter.x = ni; newinter.y = nj; newinter.z = nk;
    mListNewInter.push_back(newinter);
}

void LbmFsgrSolver::reinitFlags( int workSet ) { 
    // reinitCellFlags OLD mods:
    // add all to intel list?
    // check ffrac for new cells
    // new if cell inits (last loop)
    // vweights handling

    const int debugFlagreinit = 0;
    
    // some things need to be read/modified on the other set
    int otherSet = (workSet^1);
    // fixed level on which to perform 
    int workLev = mMaxRefine;

  /* modify interface cells from lists */
  /* mark filled interface cells as fluid, or emptied as empty */
    /* count neighbors and distribute excess mass to interface neighbor cells
   * problems arise when there are no interface neighbors anymore
     * then just distribute to any fluid neighbors...
     */

    // for symmetry, first init all neighbor cells */
    for( vector<LbmPoint>::iterator iter=mListFull.begin();
       iter != mListFull.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
        if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
    FORDF1 {
            int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
            //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
            if( (ni<=0) || (nj<=0) || 
                  (ni>=mLevel[workLev].lSizex-1) ||
                  (nj>=mLevel[workLev].lSizey-1) 
#                   if LBMDIM==3
                  || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
#                   endif // LBMDIM==1
                 ) {
                continue; } // new bc, dont treat cells on boundary NEWBC
      if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
                
                // preinit speed, get from average surrounding cells
                // interpolate from non-workset to workset, sets are handled in function

                // new and empty interface cell, dont change old flag here!
                addToNewInterList(ni,nj,nk);

                LbmFloat avgrho = 0.0;
                LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
                interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);

                // careful with l's...
                FORDF0M { 
                    QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
                    //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
                }
                //errMsg("FNEW", PRINT_VEC(ni,nj,nk)<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<<avgrho<<" vel"<<PRINT_VEC(avgux,avguy,avguz) ); // DEBUG SYMM
                QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
                QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
                //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
                changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
                if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
      }
            /* prevent surrounding interface cells from getting removed as empty cells 
             * (also cells that are not newly inited) */
      if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
                //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
                changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
                // also add to list...
                addToNewInterList(ni,nj,nk);
            } // NEW?
    }

        // NEW? no extra loop...
        //RFLAG(workLev,i,j,k, workSet) = CFFluid;
        changeFlag(workLev,i,j,k, workSet,CFFluid);
    }

    /* remove empty interface cells that are not allowed to be removed anyway
     * this is important, otherwise the dreaded cell-type-flickering can occur! */
  for(int n=0; n<(int)mListEmpty.size(); n++) {
    int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
        if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
            // treat as "new inter"
            addToNewInterList(i,j,k);
            // remove entry
            if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
            if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
            mListEmpty.pop_back();
            n--; 
        }
    } 


    /* problems arise when adjacent cells empty&fill ->
         let fill cells+surrounding interface cells have the higher importance */
  for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
       iter != mListEmpty.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
        if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
        if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );

        /* set surrounding fluid cells to interface cells */
    FORDF1 {
            int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
      if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
                // init fluid->interface 
                //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
                changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
                /* new mass = current density */
                LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
            for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
                QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
                QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 

                // store point
                addToNewInterList(ni,nj,nk);
      }
      if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
                // test, also add to list...
                addToNewInterList(ni,nj,nk);
            } // NEW?
    }

        /* for symmetry, set our flag right now */
        changeFlag(workLev,i,j,k, workSet, CFEmpty);
        // mark cell not be changed mass... - not necessary, not in list anymore anyway!
    } // emptylist


    
    // precompute weights to get rid of order dependancies
    vector<lbmFloatSet> vWeights;
    vWeights.resize( mListFull.size() + mListEmpty.size() );
    int weightIndex = 0;
  int nbCount = 0;
    LbmFloat nbWeights[LBM_DFNUM];
    LbmFloat nbTotWeights = 0.0;
    for( vector<LbmPoint>::iterator iter=mListFull.begin();
       iter != mListFull.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
    nbCount = 0; nbTotWeights = 0.0;
    FORDF1 {
            int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
      if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
                nbCount++;
                if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
                else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
                nbTotWeights += nbWeights[l];
      } else {
                nbWeights[l] = -100.0; // DEBUG;
            }
    }
        if(nbCount>0) { 
            //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
        vWeights[weightIndex].val[0] = nbTotWeights;
        FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
        vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
        } else { 
        vWeights[weightIndex].numNbs = 0.0;
        }
        weightIndex++;
    }
  for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
       iter != mListEmpty.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
    nbCount = 0; nbTotWeights = 0.0;
    FORDF1 {
            int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
      if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
                nbCount++;
                if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
                else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
                nbTotWeights += nbWeights[l];
      } else {
                nbWeights[l] = -100.0; // DEBUG;
            }
    }
        if(nbCount>0) { 
            //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
        vWeights[weightIndex].val[0] = nbTotWeights;
        FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
        vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
        } else { 
        vWeights[weightIndex].numNbs = 0.0;
        }
        weightIndex++;
    } 
    weightIndex = 0;
    

    /* process full list entries, filled cells are done after this loop */
    for( vector<LbmPoint>::iterator iter=mListFull.begin();
       iter != mListFull.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;

        LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
    FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho

    LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
        if(vWeights[weightIndex].numNbs>0.0) {
            const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
            //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
            FORDF1 {
                int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
        if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
                    LbmFloat change = -1.0;
                    if(nbTotWeightsp>0.0) {
                        //change = massChange * ( nbWeights[l]/nbTotWeightsp );
                        change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
                    } else {
                        change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
                    }
                    QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
                }
            }
            massChange = 0.0;
        } else {
            // Problem! no interface neighbors
            mFixMass += massChange;
            //TTT mNumProblems++;
            //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass);
        }
        weightIndex++;

    // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
    QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
    QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
  } // fulllist


    /* now, finally handle the empty cells - order is important, has to be after
     * full cell handling */
  for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
       iter != mListEmpty.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;

    LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
        if(vWeights[weightIndex].numNbs>0.0) {
            const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
            //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
            FORDF1 {
                int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
        if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
                    LbmFloat change = -1.0;
                    if(nbTotWeightsp>0.0) {
                        change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
                    } else {
                        change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
                    }
                    QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
                }
            }
            massChange = 0.0;
        } else {
            // Problem! no interface neighbors
            mFixMass += massChange;
            //TTT mNumProblems++;
            //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass);
        }
        weightIndex++;
        
        // finally... make it empty 
    // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
    QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
    QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
    }
  for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
       iter != mListEmpty.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
    changeFlag(workLev,i,j,k, otherSet, CFEmpty);
    } 


    // check if some of the new interface cells can be removed again 
    // never happens !!! not necessary
    // calculate ffrac for new IF cells NEW

    // how many are really new interface cells?
    int numNewIf = 0;
  for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
       iter != mListNewInter.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
        if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
            continue; 
            // FIXME remove from list?
        }
        numNewIf++;
    }

    // redistribute mass, reinit flags
    if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
    float newIfFac = 1.0/(LbmFloat)numNewIf;
  for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
       iter != mListNewInter.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
        if((i<=0) || (j<=0) || 
             (i>=mLevel[workLev].lSizex-1) ||
             (j>=mLevel[workLev].lSizey-1) ||
             ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
             ) {
            continue; } // new bc, dont treat cells on boundary NEWBC
        if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
            //errMsg("???"," "<<PRINT_IJK);
            continue; 
        } // */

    QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac);
        
        int nbored = 0;
        FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
        if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
        if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
        if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }

        if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
            RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
        }
        if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
    }

    // reinit fill fraction
  for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
       iter != mListNewInter.end(); iter++ ) {
    int i=iter->x, j=iter->y, k=iter->z;
        if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }

        initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) {
        //LbmFloat nrho = 0.0;
        //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
    //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
    //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
    }

    if(mListNewInter.size()>0){ 
        //errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() );
        mFixMass = 0.0; 
    }

    // empty lists for next step
    mListFull.clear();
    mListEmpty.clear();
    mListNewInter.clear();
} // reinitFlags