solver_init.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"
// for geo init FGI_ defines
#include "elbeem.h"

// helper for 2d init
#define SWAPYZ(vec) { \
                const LbmFloat tmp = (vec)[2]; \
                (vec)[2] = (vec)[1]; (vec)[1] = tmp; }


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

/*****************************************************************************/
#if LBMDIM==3

        const int LbmFsgrSolver::cDimension = 3;

        // Wi factors for collide step 
        const LbmFloat LbmFsgrSolver::cCollenZero    = (1.0/3.0);
        const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/18.0);
        const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);

        const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
        const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
        const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
        const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;

        const int    LbmFsgrSolver::cDfNum      = 19;
        const int    LbmFsgrSolver::cDirNum     = 27;

        //const string LbmFsgrSolver::dfString[ cDfNum ] = { 
        const char* LbmFsgrSolver::dfString[ cDfNum ] = { 
                " C", " N"," S"," E"," W"," T"," B",
                "NE","NW","SE","SW",
                "NT","NB","ST","SB",
                "ET","EB","WT","WB"
        };

        const int LbmFsgrSolver::dfNorm[ cDfNum ] = { 
                cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB, 
                cDirNE, cDirNW, cDirSE, cDirSW, 
                cDirNT, cDirNB, cDirST, cDirSB, 
                cDirET, cDirEB, cDirWT, cDirWB
        };

        const int LbmFsgrSolver::dfInv[ cDfNum ] = { 
                cDirC,  cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
                cDirSW, cDirSE, cDirNW, cDirNE,
                cDirSB, cDirST, cDirNB, cDirNT, 
                cDirWB, cDirWT, cDirEB, cDirET
        };

        const int LbmFsgrSolver::dfRefX[ cDfNum ] = { 
                0,  0, 0, 0, 0, 0, 0,
                cDirSE, cDirSW, cDirNE, cDirNW,
                0, 0, 0, 0, 
                cDirEB, cDirET, cDirWB, cDirWT
        };

        const int LbmFsgrSolver::dfRefY[ cDfNum ] = { 
                0,  0, 0, 0, 0, 0, 0,
                cDirNW, cDirNE, cDirSW, cDirSE,
                cDirNB, cDirNT, cDirSB, cDirST,
                0, 0, 0, 0
        };

        const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { 
                0,  0, 0, 0, 0, 0, 0,
                0, 0, 0, 0, 
                cDirST, cDirSB, cDirNT, cDirNB,
                cDirWT, cDirWB, cDirET, cDirEB
        };

        // Vector Order 3D:
        //  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
        //  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
        //  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
        //  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1

        const int LbmFsgrSolver::dfVecX[ cDirNum ] = { 
                0, 0,0, 1,-1, 0,0,
                1,-1,1,-1,
                0,0,0,0,
                1,1,-1,-1,
                 1,-1, 1,-1,
                 1,-1, 1,-1,
        };
        const int LbmFsgrSolver::dfVecY[ cDirNum ] = { 
                0, 1,-1, 0,0,0,0,
                1,1,-1,-1,
                1,1,-1,-1,
                0,0,0,0,
                 1, 1,-1,-1,
                 1, 1,-1,-1
        };
        const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { 
                0, 0,0,0,0,1,-1,
                0,0,0,0,
                1,-1,1,-1,
                1,-1,1,-1,
                 1, 1, 1, 1,
                -1,-1,-1,-1
        };

        const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
                0, 0,0, 1,-1, 0,0,
                1,-1,1,-1,
                0,0,0,0,
                1,1,-1,-1,
                 1,-1, 1,-1,
                 1,-1, 1,-1
        };
        const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
                0, 1,-1, 0,0,0,0,
                1,1,-1,-1,
                1,1,-1,-1,
                0,0,0,0,
                 1, 1,-1,-1,
                 1, 1,-1,-1
        };
        const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
                0, 0,0,0,0,1,-1,
                0,0,0,0,
                1,-1,1,-1,
                1,-1,1,-1,
                 1, 1, 1, 1,
                -1,-1,-1,-1
        };

        /* principal directions */
        const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { 
                1,-1, 0,0, 0,0
        };
        const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { 
                0,0, 1,-1, 0,0
        };
        const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { 
                0,0, 0,0, 1,-1
        };

        LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
        LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];


        const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { 
                cCollenZero,
                cCollenOne, cCollenOne, cCollenOne, 
                cCollenOne, cCollenOne, cCollenOne,
                cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
                cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
                cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
        };

        /* precalculated equilibrium dfs, inited in lbmsolver constructor */
        LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];

#else // end LBMDIM==3 , LBMDIM==2

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

                const int LbmFsgrSolver::cDimension = 2;

                const LbmFloat LbmFsgrSolver::cCollenZero    = (4.0/9.0);
                const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/9.0);
                const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);

                const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
                const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
                const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
                const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;

                const int LbmFsgrSolver::cDfNum  = 9;
                const int LbmFsgrSolver::cDirNum = 9;

        //const string LbmFsgrSolver::dfString[ cDfNum ] = { 
        const char* LbmFsgrSolver::dfString[ cDfNum ] = { 
                " C", 
                " N",   " S", " E", " W",
                "NE", "NW", "SE","SW" 
        };

        const int LbmFsgrSolver::dfNorm[ cDfNum ] = { 
                cDirC, 
                cDirN,  cDirS,  cDirE,  cDirW,
                cDirNE, cDirNW, cDirSE, cDirSW 
        };

        const int LbmFsgrSolver::dfInv[ cDfNum ] = { 
                cDirC,  
                cDirS,  cDirN,  cDirW,  cDirE,
                cDirSW, cDirSE, cDirNW, cDirNE 
        };

        const int LbmFsgrSolver::dfRefX[ cDfNum ] = { 
                0,  
                0,  0,  0,  0,
                cDirSE, cDirSW, cDirNE, cDirNW 
        };

        const int LbmFsgrSolver::dfRefY[ cDfNum ] = { 
                0,  
                0,  0,  0,  0,
                cDirNW, cDirNE, cDirSW, cDirSE 
        };

        const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { 
                0,  0, 0, 0, 0,
                0, 0, 0, 0
        };

        // Vector Order 2D:
        // 0  1 2  3  4  5  6 7  8
        // 0, 0,0, 1,-1, 1,-1,1,-1 
        // 0, 1,-1, 0,0, 1,1,-1,-1 

        const int LbmFsgrSolver::dfVecX[ cDirNum ] = { 
                0, 
                0,0, 1,-1,
                1,-1,1,-1 
        };
        const int LbmFsgrSolver::dfVecY[ cDirNum ] = { 
                0, 
                1,-1, 0,0,
                1,1,-1,-1 
        };
        const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { 
                0, 0,0,0,0, 0,0,0,0 
        };

        const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
                0, 
                0,0, 1,-1,
                1,-1,1,-1 
        };
        const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
                0, 
                1,-1, 0,0,
                1,1,-1,-1 
        };
        const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
                0, 0,0,0,0, 0,0,0,0 
        };

        const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { 
                1,-1, 0,0
        };
        const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { 
                0,0, 1,-1
        };
        const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { 
                0,0, 0,0
        };


        LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
        LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];


        const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { 
                cCollenZero,
                cCollenOne, cCollenOne, cCollenOne, cCollenOne, 
                cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
        };

        /* precalculated equilibrium dfs, inited in lbmsolver constructor */
        LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];

// D2Q9 end
#endif  // LBMDIM==2


// required globals
extern bool glob_mpactive;
extern int glob_mpnum, glob_mpindex;


/******************************************************************************
 * Lbm Constructor
 *****************************************************************************/
LbmFsgrSolver::LbmFsgrSolver() :
        //D(),
        mCurrentMass(0.0), mCurrentVolume(0.0),
        mNumProblems(0), 
        mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
        mpPreviewSurface(NULL), 
        mTimeAdap(true), mForceTimeStepReduce(false),
        mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
        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),
        mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
        mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
        mIsoWeightMethod(1),
        mMaxRefine(1), 
        mDfScaleUp(-1.0), mDfScaleDown(-1.0),
        mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03
        mDebugOmegaRet(0.0),
        mLastOmega(1e10), mLastGravity(1e10),
        mNumInvIfTotal(0), mNumFsgrChanges(0),
        mDisableStandingFluidInit(0),
        mInit2dYZ(false),
        mForceTadapRefine(-1), mCutoff(-1)
{
        mpControl = new LbmControlData();

#if LBM_INCLUDE_TESTSOLVERS==1
        mpTest = new LbmTestdata();
        mMpNum = mMpIndex = 0;
        mOrgSizeX = 0;
        mOrgStartX = 0.;
        mOrgEndX = 0.;
#endif // LBM_INCLUDE_TESTSOLVERS!=1
        mpIso = new IsoSurface( mIsoValue );

  // init equilibrium dist. func 
  LbmFloat rho=1.0;
  FORDF0 {
                dfEquil[l] = this->getCollideEq( l,rho,  0.0, 0.0, 0.0);
  }
        dfEquil[dMass] = 1.;
        dfEquil[dFfrac] = 1.;
        dfEquil[dFlux] = FLUX_INIT;

        // init LES
        int odm = 0;
        for(int m=0; m<LBMDIM; m++) { 
                for(int l=0; l<cDfNum; l++) { 
                        this->lesCoeffDiag[m][l] = 
                        this->lesCoeffOffdiag[m][l] = 0.0;
                }
        }
        for(int m=0; m<LBMDIM; m++) { 
                for(int n=0; n<LBMDIM; n++) { 
                        for(int l=1; l<cDfNum; l++) { 
                                LbmFloat em;
                                switch(m) {
                                        case 0: em = dfDvecX[l]; break;
                                        case 1: em = dfDvecY[l]; break;
                                        case 2: em = dfDvecZ[l]; break;
                                        default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
                                }
                                LbmFloat en;
                                switch(n) {
                                        case 0: en = dfDvecX[l]; break;
                                        case 1: en = dfDvecY[l]; break;
                                        case 2: en = dfDvecZ[l]; break;
                                        default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
                                }
                                const LbmFloat coeff = em*en;
                                if(m==n) {
                                        this->lesCoeffDiag[m][l] = coeff;
                                } else {
                                        if(m>n) {
                                                this->lesCoeffOffdiag[odm][l] = coeff;
                                        }
                                }
                        }

                        if(m==n) {
                        } else {
                                if(m>n) odm++;
                        }
                }
        }

        mDvecNrm[0] = LbmVec(0.0);
  FORDF1 {
                mDvecNrm[l] = getNormalized( 
                        LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] ) 
                        ) * -1.0; 
        }

        // calculate gauss weights for restriction
        //LbmFloat mGaussw[27];
        LbmFloat totGaussw = 0.0;
        const LbmFloat alpha = 1.0;
        const LbmFloat gw = sqrt(2.0*LBMDIM);
#if ELBEEM_PLUGIN!=1
        errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
#endif
        for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
        //for(int n=0;(n<cDirNum); n++) { 
        for(int n=0;(n<cDfNum); n++) { 
                const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
                LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
                mGaussw[n] = w;
                totGaussw += w;
        }
        for(int n=0;(n<cDirNum); n++) { 
                mGaussw[n] = mGaussw[n]/totGaussw;
        }

}

/*****************************************************************************/
/* Destructor */
/*****************************************************************************/
LbmFsgrSolver::~LbmFsgrSolver()
{
  if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; }
#if COMPRESSGRIDS==1
        delete [] mLevel[mMaxRefine].mprsCells[1];
        mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
#endif // COMPRESSGRIDS==1

        for(int i=0; i<=mMaxRefine; i++) {
                for(int s=0; s<2; s++) {
                        if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
                        if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
                }
        }
        delete mpIso;
        if(mpPreviewSurface) delete mpPreviewSurface;
        // cleanup done during scene deletion...
        
        if(mpControl) delete mpControl;

        // always output performance estimate
        debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
  if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
}




/******************************************************************************
 * initilize variables fom attribute list 
 *****************************************************************************/
void LbmFsgrSolver::parseAttrList()
{
        LbmSolverInterface::parseStdAttrList();

        string matIso("default");
        matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
        mpIso->setMaterialName( matIso );
        mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
        mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
        mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false );
        mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );

        mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
        mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
        mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
        mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
        mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );

        // refinement
        mMaxRefine = mRefinementDesired;
        mMaxRefine  = mpSifAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
        if(mMaxRefine<0) mMaxRefine=0;
        if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
        mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
        mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
        mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);

        // demo mode settings
        mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
        // FIXME check needed?
        mFVArea   = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );

        // debugging - skip some time...
        double starttimeskip = 0.;
        starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
        mSimulationTime += starttimeskip;
        if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);

        mpControl->parseControldataAttrList(mpSifAttrs);

#if LBM_INCLUDE_TESTSOLVERS==1
        mUseTestdata = 0;
        mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
        mpTest->parseTestdataAttrList(mpSifAttrs);
#ifdef ELBEEM_PLUGIN
        mUseTestdata=1; // DEBUG
#endif // ELBEEM_PLUGIN

        mMpNum = mpSifAttrs->readInt("mpnum",  mMpNum ,"LbmFsgrSolver", "mMpNum", false);
        mMpIndex = mpSifAttrs->readInt("mpindex",  mMpIndex ,"LbmFsgrSolver", "mMpIndex", false);
        if(glob_mpactive) {
                // used instead...
                mMpNum = glob_mpnum;
                mMpIndex = glob_mpindex;
        } else {
                glob_mpnum = mMpNum;
                glob_mpindex = 0;
        }
        errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<<glob_mpactive<<", "<<glob_mpindex<<"/"<<glob_mpnum);
        if(mMpNum>0) {
                mUseTestdata=1; // needed in this case...
        }

        errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
#else // LBM_INCLUDE_TESTSOLVERS!=1
        // not testsolvers, off by default
        mUseTestdata = 0;
        if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
#endif // LBM_INCLUDE_TESTSOLVERS!=1

        mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
        // deprecated!
        float mInitialCsmagoCoarse = 0.0;
        mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
#if USE_LES==1
#else // USE_LES==1
        debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
        mInitialCsmago = 0.0;
#endif // USE_LES==1
}


/******************************************************************************
 * (part of enabling chapter 6 of "Free Surface Flows with Moving and Deforming Objects for LBM")
 *****************************************************************************/
void LbmFsgrSolver::setSurfGenSettings(short value)
{
        mFsSurfGenSetting = value;
}


/******************************************************************************
 * Initialize omegas and forces on all levels (for init/timestep change)
 *****************************************************************************/
void LbmFsgrSolver::initLevelOmegas()
{
        // no explicit settings
        mOmega = mpParam->calculateOmega(mSimulationTime);
        mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) );
        mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused
        if(mInit2dYZ) { SWAPYZ(mGravity); }

        // check if last init was ok
        LbmFloat gravDelta = norm(mGravity-mLastGravity);
        //errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<mOmega<<" - lom:"<<mLastOmega<<" gv:"<<mGravity<<" - "<<mLastGravity<<" , "<<gravDelta  );
        if((mOmega == mLastOmega) && (gravDelta<=0.0)) return;

        if(mInitialCsmago<=0.0) {
                if(OPT3D==1) {
                        errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR); 
                        return;
                }
        }

        LbmFloat  fineCsmago  = mInitialCsmago;
        LbmFloat coarseCsmago = mInitialCsmago;
        LbmFloat maxFineCsmago1    = 0.026; 
        LbmFloat maxCoarseCsmago1  = 0.029; // try stabilizing
        LbmFloat maxFineCsmago2    = 0.028; 
        LbmFloat maxCoarseCsmago2  = 0.032; // try stabilizing some more
        if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) {
                fineCsmago = maxFineCsmago1;
                coarseCsmago = maxCoarseCsmago1;
        }
        if((mMaxRefine>1)&&(mInitialCsmago<maxFineCsmago2)) {
                fineCsmago = maxFineCsmago2;
                coarseCsmago = maxCoarseCsmago2;
        }
                

        // use Tau instead of Omega for calculations
        { // init base level
                int i = mMaxRefine;
                mLevel[i].omega    = mOmega;
                mLevel[i].timestep = mpParam->getTimestep();
                mLevel[i].lcsmago = fineCsmago; //CSMAGO_INITIAL;
                mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
                mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
        }

        // init all sub levels
        for(int i=mMaxRefine-1; i>=0; i--) {
                //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
                double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
                nomega                = 1.0/nomega;
                mLevel[i].omega       = (LbmFloat)nomega;
                mLevel[i].timestep    = 2.0 * mLevel[i+1].timestep;
                mLevel[i].lcsmago = coarseCsmago;
                mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
                mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
        }
        
        // for lbgk
        mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
        for(int i=mMaxRefine-1; i>=0; i--) {
                // should be the same on all levels...
                // for lbgk
                mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
        }

        mLastOmega = mOmega;
        mLastGravity = mGravity;
        // debug? invalidate old values...
        mGravity = -100.0;
        mOmega = -100.0;

        for(int i=0; i<=mMaxRefine; i++) {
                if(!mSilent) {
                        errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz 
                                        <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
                                        <<" cmsagp:"<<mLevel[i].lcsmago<<", "
                                        << " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
                } else {
                        if(!mInitDone) {
                                debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
                                                <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
                        }
                }
        }
        if(mMaxRefine>0) {
                mDfScaleUp   = (mLevel[0  ].timestep/mLevel[0+1].timestep)* (1.0/mLevel[0  ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
                mDfScaleDown = (mLevel[0+1].timestep/mLevel[0  ].timestep)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0  ].omega-1.0); // yu
        }
}


/******************************************************************************
 * Init Solver (values should be read from config file)
 *****************************************************************************/

bool LbmFsgrSolver::initializeSolverMemory()
{
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<mInitDone<<" "<<(void*)this,1);

        // init cppf stage
        if(mCppfStage>0) {
                mSizex *= mCppfStage;
                mSizey *= mCppfStage;
                mSizez *= mCppfStage;
        }
        if(mFsSurfGenSetting==-1) {
                // all on
                mFsSurfGenSetting = 
                         fssgNormal   | fssgNoNorth  | fssgNoSouth  | fssgNoEast   |
                         fssgNoWest   | fssgNoTop    | fssgNoBottom | fssgNoObs   ;
        }

        // size inits to force cubic cells and mult4 level dimensions
        // and make sure we dont allocate too much...
        bool memOk = false;
        int orgSx = mSizex;
        int orgSy = mSizey;
        int orgSz = mSizez;
        double sizeReduction = 1.0;
        double memEstFromFunc = -1.0;
        double memEstFine = -1.0;
        string memreqStr("");   
        bool firstMInit = true;
        int minitTries=0;
        while(!memOk) {
                minitTries++;
                initGridSizes( mSizex, mSizey, mSizez,
                                mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL);

                // MPT
#if LBM_INCLUDE_TESTSOLVERS==1
                if(firstMInit) {
                        mrSetup();
                }
#endif // LBM_INCLUDE_TESTSOLVERS==1
                firstMInit=false;

                calculateMemreqEstimate( mSizex, mSizey, mSizez, 
                                mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr );
                
                double memLimit;
                string memLimStr("-");
                if(sizeof(void*)==4) {
                        // 32bit system, limit to 2GB
                        memLimit = 2.0* 1024.0*1024.0*1024.0;
                        memLimStr = string("2GB");
                } else {
                        // 64bit, just take 16GB as limit for now...
                        memLimit = 16.0* 1024.0*1024.0*1024.0;
                        memLimStr = string("16GB");
                }

                // restrict max. chunk of 1 mem block to 1GB for windos
                bool memBlockAllocProblem = false;
                double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
                //std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
#ifdef WIN32
                double maxWinMemChunk = 1100.*1024.*1024.;
                if(sizeof(void *)==4 && memEstFine>maxWinMemChunk) {
                        memBlockAllocProblem = true;
                }
#endif // WIN32
#ifdef __APPLE__
                double maxMacMemChunk = 1200.*1024.*1024.;
                if(memEstFine> maxMacMemChunk) {
                        memBlockAllocProblem = true;
                }
#endif // Mac
                if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
                        // max memory chunk for 32bit systems 2gig
                        memBlockAllocProblem = true;
                }

                if(memEstFromFunc>memLimit || memBlockAllocProblem) {
                        sizeReduction *= 0.9;
                        mSizex = (int)(orgSx * sizeReduction);
                        mSizey = (int)(orgSy * sizeReduction);
                        mSizez = (int)(orgSz * sizeReduction);
                        debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
                                        //memEstFromFunc<<"/"<<memLimit<<", "<<
                                        //memEstFine<<"/"<<maxWinMemChunk<<", "<<
                                        memreqStr<<"/"<<memLimStr<<", "<<
                                        "retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
                                        , 3 );
                } else {
                        memOk = true;
                } 
        }
        
        mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
          ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
          ", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
          ", est. Mem.Req.: "<<memreqStr        ,2);
        mpParam->setSize(mSizex, mSizey, mSizez);
        if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }

  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
                <<"LBM_EPSILON="<<LBM_EPSILON  <<" "
                <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG <<" "
                <<"OPT3D="<<OPT3D <<" "
                <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
                <<"MASS_INVALID="<<MASS_INVALID <<" "
                <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
                <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
                <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
                <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" " 
                <<"USE_LES="<<USE_LES <<" " 
                ,10);

        // perform 2D corrections...
        if(LBMDIM == 2) mSizez = 1;

        mpParam->setSimulationMaxSpeed(0.0);
        if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
        mpParam->setTadapLevels( mMaxRefine+1 );

        if(mForceTadapRefine>mMaxRefine) {
                mpParam->setTadapLevels( mForceTadapRefine+1 );
                debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
        }

        if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
                errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
                return false;
        }


        // init vectors
        for(int i=0; i<=mMaxRefine; i++) {
                mLevel[i].id = i;
                mLevel[i].nodeSize = 0.0; 
                mLevel[i].simCellSize = 0.0; 
                mLevel[i].omega = 0.0; 
                mLevel[i].time = 0.0; 
                mLevel[i].timestep = 1.0;
                mLevel[i].gravity = LbmVec(0.0); 
                mLevel[i].mprsCells[0] = NULL;
                mLevel[i].mprsCells[1] = NULL;
                mLevel[i].mprsFlags[0] = NULL;
                mLevel[i].mprsFlags[1] = NULL;

                mLevel[i].avgOmega = 0.0; 
                mLevel[i].avgOmegaCnt = 0.0;
        }

        // init sizes
        mLevel[mMaxRefine].lSizex = mSizex;
        mLevel[mMaxRefine].lSizey = mSizey;
        mLevel[mMaxRefine].lSizez = mSizez;
        for(int i=mMaxRefine-1; i>=0; i--) {
                mLevel[i].lSizex = mLevel[i+1].lSizex/2;
                mLevel[i].lSizey = mLevel[i+1].lSizey/2;
                mLevel[i].lSizez = mLevel[i+1].lSizez/2;
        }

        // safety check
        if(sizeof(CellFlagType) != CellFlagTypeSize) {
                errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
                return false;
        }

        double ownMemCheck = 0.0;
        mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
        mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
        mLevel[ mMaxRefine ].lcellfactor = 1.0;
        LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);

#if COMPRESSGRIDS==0
        mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
        mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
        ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
#else // COMPRESSGRIDS==0
        LONGINT compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
        // D int tmp = ( (rcellSize +compressOffset +4)/(1024*1024) )*4;
        // D printf("Debug MEMMMM excee: %d\n", tmp);
        mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
        mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
        ownMemCheck += sizeof(LbmFloat) * (rcellSize +compressOffset +4);
#endif // COMPRESSGRIDS==0

        if(!mLevel[ mMaxRefine ].mprsCells[1] || !mLevel[ mMaxRefine ].mprsCells[0]) {
                errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (1)! Aborting...",SIMWORLD_INITERROR);
                return false;
        }

        // +4 for safety ?
        mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
        mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
        ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
        if(!mLevel[ mMaxRefine ].mprsFlags[1] || !mLevel[ mMaxRefine ].mprsFlags[0]) {
                errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (2)! Aborting...",SIMWORLD_INITERROR);

#if COMPRESSGRIDS==0
                delete[] mLevel[ mMaxRefine ].mprsCells[0];
                delete[] mLevel[ mMaxRefine ].mprsCells[1];
#else // COMPRESSGRIDS==0
                delete[] mLevel[ mMaxRefine ].mprsCells[1];
#endif // COMPRESSGRIDS==0
                return false;
        }

        LbmFloat lcfdimFac = 8.0;
        if(LBMDIM==2) lcfdimFac = 4.0;
        for(int i=mMaxRefine-1; i>=0; i--) {
                mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
                mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
                mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;

                if(LBMDIM==2){ mLevel[i].lSizez = 1; } // 2D
                rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
                mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
                mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
                ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
                mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
                mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
                ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
        }

        // isosurface memory, use orig res values
        if(mFarFieldSize>0.) {
                ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
        } else {
                // ignore 3 int slices...
                ownMemCheck += (double)( (              sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
        }

        // sanity check
#if ELBEEM_PLUGIN!=1
        if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
                errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
        }
#endif // ELBEEM_PLUGIN!=1
        
        // init sizes for _all_ levels
        for(int i=mMaxRefine; i>=0; i--) {
                mLevel[i].lOffsx = mLevel[i].lSizex;
                mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
                mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
        mLevel[i].setCurr  = 0;
        mLevel[i].setOther = 1;
        mLevel[i].lsteps = 0;
        mLevel[i].lmass = 0.0;
        mLevel[i].lvolume = 0.0;
        }

        // calc omega, force for all levels
        initLevelOmegas();
        mMinTimestep = mpParam->getTimestep();
        mMaxTimestep = mpParam->getTimestep();

        // init isosurf
        mpIso->setIsolevel( mIsoValue );
#if LBM_INCLUDE_TESTSOLVERS==1
        if(mUseTestdata) {
                mpTest->setMaterialName( mpIso->getMaterialName() );
                delete mpIso;
                mpIso = mpTest;
                if(mpTest->mFarfMode>0) { // 3d off
                        mpTest->setIsolevel(-100.0);
                } else {
                        mpTest->setIsolevel( mIsoValue );
                }
        }
#endif // LBM_INCLUDE_TESTSOLVERS!=1
        // approximate feature size with mesh resolution
        float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
        // smooth vars defined in solver_interface, set by simulation object
        // reset for invalid values...
        if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
        if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
        mpIso->setSmoothSurface( mSmoothSurface * featureSize );
        mpIso->setSmoothNormals( mSmoothNormals * featureSize );

        // init iso weight values mIsoWeightMethod
        int wcnt = 0;
        float totw = 0.0;
        for(int ak=-1;ak<=1;ak++) 
                for(int aj=-1;aj<=1;aj++) 
                        for(int ai=-1;ai<=1;ai++)  {
                                switch(mIsoWeightMethod) {
                                case 1: // light smoothing
                                        mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
                                        break;
                                case 2: // very light smoothing
                                        mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
                                        mIsoWeight[wcnt] *= mIsoWeight[wcnt];
                                        break;
                                case 3: // no smoothing
                                        if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
                                        else mIsoWeight[wcnt] = 0.0;
                                        break;
                                default: // strong smoothing (=0)
                                        mIsoWeight[wcnt] = 1.0;
                                        break;
                                }
                                totw += mIsoWeight[wcnt];
                                wcnt++;
                        }
        for(int i=0; i<27; i++) mIsoWeight[i] /= totw;

        LbmVec isostart = vec2L(mvGeoStart);
        LbmVec isoend   = vec2L(mvGeoEnd);
        int twodOff = 0; // 2d slices
        if(LBMDIM==2) {
                LbmFloat sn,se;
                sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
                se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
                isostart[2] = sn;
                isoend[2] = se;
                twodOff = 2;
        }
        int isosx = mSizex+2;
        int isosy = mSizey+2;
        int isosz = mSizez+2+twodOff;

        // MPT
#if LBM_INCLUDE_TESTSOLVERS==1
        //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
        if( (mMpNum>0) && (mMpIndex==0) ) {
                //? mpindex==0
                // restore original value for node0
                isosx       = mOrgSizeX + 2;
                isostart[0] = mOrgStartX;
                isoend[0]   = mOrgEndX;
        }
        errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
                        << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
                        << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
        errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
                        << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
                        << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
#endif // LBM_INCLUDE_TESTSOLVERS==1

        errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
        errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
                        " grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
        mpIso->setStart( vec2G(isostart) );
        mpIso->setEnd(   vec2G(isoend) );
        LbmVec isodist = isoend-isostart;

        int isosubs = mIsoSubdivs;
        if(mFarFieldSize>1.) {
                errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
                isosubs = 1;
                mpIso->setUseFulledgeArrays(true);
        }
        mpIso->setSubdivs(isosubs);

        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) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
                        if(!mAllfluid) {
                                initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
                        } else {
                                initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); 
                        }
                }
        }


        if(LBMDIM==2) {
                if(mOutputSurfacePreview) {
                        errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
                        mOutputSurfacePreview = 0; }
        }
        if((glob_mpactive) && (glob_mpindex>0)) {
                mOutputSurfacePreview = 0;
        }

#if LBM_USE_GUI==1
        if(mOutputSurfacePreview) {
                errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
                mOutputSurfacePreview = 0; }
#endif // LBM_USE_GUI==1
        if(mOutputSurfacePreview) {
                // same as normal one, but use reduced size
                mpPreviewSurface = new IsoSurface( mIsoValue );
                mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
                mpPreviewSurface->setIsolevel( mIsoValue );
                // usually dont display for rendering
                mpPreviewSurface->setVisible( false );

                mpPreviewSurface->setStart( vec2G(isostart) );
                mpPreviewSurface->setEnd(   vec2G(isoend) );
                LbmVec pisodist = isoend-isostart;
                LbmFloat pfac = mPreviewFactor;
                mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
                //mpPreviewSurface->setName( getName() + "preview" );
                mpPreviewSurface->setName( "preview" );
        
                debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
        }

        // init defaults
        mAvgNumUsedCells = 0;
        mFixMass= 0.0;
        return true;
}

bool LbmFsgrSolver::initializeSolverGrids() {
  /* init boundaries */
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
        // init obstacles, and reinit time step size 
        initGeometryFlags();
        mLastSimTime = -1.0;
        // TODO check for invalid cells? nitGenericTestCases();
        
        // new - init noslip 1 everywhere...
        // half fill boundary cells?

        CellFlagType domainBoundType = CFInvalid;
        // TODO use normal object types instad...
        if(mDomainBound.find(string("free")) != string::npos) {
                domainBoundType = CFBnd | CFBndFreeslip;
        debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
        } else if(mDomainBound.find(string("part")) != string::npos) {
                domainBoundType = CFBnd | CFBndPartslip; // part slip type
        debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
        } else { 
                domainBoundType = CFBnd | CFBndNoslip;
        debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
        }

        // use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
        int domainobj = (int)(mpGiObjects->size());
        domainBoundType |= (domainobj<<24);
        //for(int i=0; i<(int)(domainobj+0); i++) {
                //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
                //if((*mpGiObjects)[i] == mpIso) { //check...
                //} 
        //}
        //errMsg("GEOIN"," dm "<<(domainBoundType>>24));

  for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
    for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
                        initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL); 
                        initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL); 
    }

  for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
    for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
                        initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL); 
                        initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL); 
                        // DEBUG BORDER!
                        //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
    }

        if(LBMDIM == 3) {
                // only for 3D
                for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
                        for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
                                initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL); 
                                initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL); 
                        }
        }

        // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
  /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
    for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
                        initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
    }
  for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
    for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
                        initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL); 
    }
        // */

        /*for(int ii=0; ii<(int)po w_change?(2.0,mMaxRefine)-1; ii++) {
                errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
                for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
                        for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
                                initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL);  // SYMM!? 2D?
                        }
                for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
                        for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
                                initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL);   // SYMM!? 3D?
                        }
        }
        // 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; }

        // prepare interface cells
        initFreeSurfaces();
        initStandingFluidGradient();

        // perform first step to init initial mass
        mInitialMass = 0.0;
        int inmCellCnt = 0;
        FSGR_FORIJK1(mMaxRefine) {
                if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
                        LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); 
                        FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
                        mInitialMass += fluidRho;
                        inmCellCnt ++;
                } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
                        mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
                        inmCellCnt ++;
                }
        }
        mCurrentVolume = mCurrentMass = mInitialMass;

        ParamVec cspv = mpParam->calculateCellSize();
        if(LBMDIM==2) cspv[2] = 1.0;
        inmCellCnt = 1;
        double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
        debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
        mInitialMass = 0.0; // reset, and use actual value after first step

        //mStartSymm = false;
#if ELBEEM_PLUGIN!=1
        if((LBMDIM==2)&&(mSizex<200)) {
                if(!checkSymmetry("init")) {
                        errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
                } else {
                        errMsg("LbmFsgrSolver::initialize","Symmetric init!");
                }
        }
#endif // ELBEEM_PLUGIN!=1
        return true;
}


bool LbmFsgrSolver::initializeSolverPostinit() {
        // coarsen region
        myTime_t fsgrtstart = getTime(); 
        for(int lev=mMaxRefine-1; lev>=0; lev--) {
                debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
                adaptGrid(lev);
                coarseRestrictFromFine(lev);
                adaptGrid(lev);
                coarseRestrictFromFine(lev);
        }
        markedClearList();
        myTime_t fsgrtend = getTime(); 
        if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
        mNumFsgrChanges = 0;

        for(int l=0; l<cDirNum; l++) { 
                LbmFloat area = 0.5 * 0.5 *0.5;
                if(LBMDIM==2) area = 0.5 * 0.5;

                if(dfVecX[l]!=0) area *= 0.5;
                if(dfVecY[l]!=0) area *= 0.5;
                if(dfVecZ[l]!=0) area *= 0.5;
                mFsgrCellArea[l] = area;
        } // l

        // make sure both sets are ok
        // copy from other to curr
        for(int lev=0; lev<=mMaxRefine; lev++) {
        FSGR_FORIJK_BOUNDS(lev) {
                RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
        } } // first copy flags */


        // old mpPreviewSurface init
        //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
        // make sure fill fracs are right for first surface generation
        stepMain();

        // prepare once...
        mpIso->setParticles(mpParticles, mPartDropMassSub);
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
        prepareVisualization();
        // copy again for stats counting
        for(int lev=0; lev<=mMaxRefine; lev++) {
        FSGR_FORIJK_BOUNDS(lev) {
                RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
        } } // first copy flags */


        // now really done...
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
  debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
        mInitDone = 1;

        // init fluid control
        initCpdata();

#if LBM_INCLUDE_TESTSOLVERS==1
        initTestdata();
#endif // ELBEEM_PLUGIN!=1
        // not inited? dont use...
        if(mCutoff<0) mCutoff=0;

        initParticles();
        return true;
}



// macros for mov obj init
#if LBMDIM==2

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

#else // LBMDIM -> 3
#define POS2GRID_CHECK(vec,n) \
                                monTotal++;\
                                const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
                                if(i<=0) continue; \
                                if(i>=mLevel[level].lSizex-1) continue; \
                                const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
                                if(j<=0) continue; \
                                if(j>=mLevel[level].lSizey-1) continue; \
                                const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
                                if(k<=0) continue; \
                                if(k>=mLevel[level].lSizez-1) continue; \

#endif // LBMDIM 

// calculate object velocity from vert arrays in objvel vec
#define OBJVEL_CALC \
                                LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
                                const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
                                USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
                                if(usqr>maxusqr) { \
                                        /* cutoff at maxVelVal */ \
                                        for(int jj=0; jj<3; jj++) { \
                                                if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
                                                if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
                                        } \
                                } } \
                                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, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
                                        /* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
                                        objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
                                }

#define TTT \


/*****************************************************************************/
/*****************************************************************************/
void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
        myTime_t monstart = getTime();

        // 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 + mpParam->getTimestep();
        if(mLastSimTime == targetTime) {
                debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
                return;
        }
        //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
        //if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);

        const LbmFloat maxVelVal = 0.1666;
        const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;

        LbmFloat rhomass = 0.0;
        CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
        CellFlagType ntype = CFInvalid;
        // WARNING - copied from geo init!
        int numobjs = (int)(mpGiObjects->size());
        ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
        // 2d display as rectangles
        ntlVec3Gfx iniPos(0.0);
        if(LBMDIM==2) {
                dvec[2] = 1.0; 
                iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
        } else {
                iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
        }
        
        if( (int)mObjectMassMovnd.size() < numobjs) {
                for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
                        mObjectMassMovnd.push_back(0.);
                }
        }
        
        // stats
        int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
        int nbored;
        for(int OId=0; OId<numobjs; OId++) {
                ntlGeometryObject *obj = (*mpGiObjects)[OId];
                bool skip = false;
                if(obj->getGeoInitId() != mLbmInitId) skip=true;
                if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
                if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
                if(skip) continue;
                debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);

                if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || 
                                ((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) {

                        otype = ntype = CFInvalid;
                        switch(obj->getGeoInitType()) {
                                /* case FGI_BNDPART: // old, use noslip for moving part/free objs
                                case FGI_BNDFREE: 
                                        if(!staticInit) {
                                                errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
                                                otype = ntype = CFBnd|CFBndNoslip;
                                        } else {
                                                if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
                                                if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
                                        }
                                        break; 
                                        // off */
                                case FGI_BNDPART: rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
                                        break;
                                case FGI_BNDFREE: rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
                                        break;
                                        // off */
                                case FGI_BNDNO:   rhomass = BND_FILL;
                                        otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
                                        break;
                                case FGI_FLUID: 
                                        otype = ntype = CFFluid; 
                                        break;
                                case FGI_MBNDINFLOW: 
                                        otype = ntype = CFMbndInflow; 
                                        break;
                                case FGI_MBNDOUTFLOW: 
                                        otype = ntype = CFMbndOutflow; 
                                        break;
                        }
                        int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
                        int active =    ((obj->getGeoActive(targetTime)>0.)? 1:0);
                        //errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<"  s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
                        // skip inactive in/out flows
                        if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
                        if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;

                        // copied from  recalculateObjectSpeeds
                        mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
                        debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*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, pNormals,  mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
                                monTrafo += mMOIVerticesOld.size();
                                obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
                                monTrafo += mMOIVertices.size();
                                obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
                        } else {
                                // only do transform update
                                obj->getMovingPoints(mMOIVertices,pNormals); // mMOIVertices = mCachedMovPoints
                                mMOIVerticesOld = mMOIVertices;
                                // WARNING - assumes mSimulationTime is global!?
                                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 /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
                                monTrafo += mMOIVerticesOld.size();
                        }

                        // object types
                        if(ntype&CFBnd){

                                // check if object is moving at all
                                if(obj->getIsAnimated()) {
                                        ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
                                        // FIXME?
                                        if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
                                        // get old type - CHECK FIXME , timestep could have changed - cause trouble?
                                        ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
                                        if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
                                }
                                if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
                                CellFlagType rflagnb[27];
                                LbmFloat massCheck = 0.;
                                int massReinits=0;                              
                                bool fillCells = (mObjectMassMovnd[OId]<=-1.);
                                LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
                                

                                // first pass - set new obs. cells
                                if(active) {
                                        for(size_t n=0; n<mMOIVertices.size(); n++) {
                                                //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
                                                POS2GRID_CHECK(mMOIVertices,n);
                                                //{ 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); }
                                                        massReinits++;
                                                }
                                                else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
                                                        massCheck -= QCELL(level, i,j,k, workSet, dMass);
                                                        massReinits++;
                                                }

                                                RFLAG(level, i,j,k, workSet) = ntype;
                                                FORDF1 {
                                                        //CellFlagType flag = RFLAG_NB(level, i,j,k,workSet,l);
                                                        rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
                                                        if(rflagnb[l]&(CFFluid|CFInter)) {
                                                                rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
                                                                RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l]; 
                                                        }
                                                }
                                                LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
                                                RAC(dstCell,0) = 0.0;
                                                if(ntype&CFBndMoving) {
                                                        OBJVEL_CALC;
                                                        
                                                        // compute fluid acceleration
                                                        FORDF1 {
                                                                if(rflagnb[l]&(CFFluid|CFInter)) { 
                                                                        const LbmFloat ux = dfDvecX[l]*objvel[0];
                                                                        const LbmFloat uy = dfDvecY[l]*objvel[1];
                                                                        const LbmFloat uz = dfDvecZ[l]*objvel[2];

                                                                        LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // 
                                                                        if(ntype&(CFBndFreeslip|CFBndPartslip)) {
                                                                                // missing, diag mass checks...
                                                                                //if(l<=LBMDIM*2) factor *= 1.4142;
                                                                                factor *= 2.0; // TODO, optimize
                                                                        } else {
                                                                                factor *= 1.2; // TODO, optimize
                                                                        }
                                                                        factor *= impactCorrFactor; // use manual tweaking channel
                                                                        RAC(dstCell,l) = factor;
                                                                        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; }
                                                }
                                                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)) {
                                                        // 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,otherSet,l); // test use other set to not have loop dependance
                                                                nbored |= rflagnb[l];
                                                        } 
                                                        CellFlagType settype = CFInvalid;
                                                        if(nbored&CFFluid) {
                                                                settype = CFInter|CFNoInterpolSrc; 
                                                                rhomass = 1.5;
                                                                if(!fillCells) rhomass = 0.;

                                                                OBJVEL_CALC;
                                                                if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }

                                                                // 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);
                                                                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) );
                                                                massCheck += rhomass;
                                                        } else {
                                                                settype = CFEmpty; rhomass = 0.0;
                                                                initEmptyCell(level, i,j,k, settype, 1., rhomass );
                                                        }
                                                        monFluids++;
                                                        massReinits++;
                                                } // flag & simtime
                                        }
                                }  // wasactive

                                // only compute mass transfer when init is done
                                if(mInitDone) {
                                        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

                        else if(ntype&CFFluid){
                                // second static init pass
                                if(staticInit) {
                                        //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" verts"<<mMOIVertices.size() ,9);
                                        CellFlagType setflflag = CFFluid; //|(OId<<24);
                                        for(size_t n=0; n<mMOIVertices.size(); n++) {
                                                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, setflflag);
                                                        initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
                                                } 
                                                //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
                                        }
                                } // second static inflow pass
                        } // fluid

                        else if(ntype&CFMbndInflow){
                                // inflow pass - add new fluid cells
                                // this is slightly different for standing inflows,
                                // as the fluid is forced to the desired velocity inside as 
                                // well...
                                const LbmFloat iniRho = 1.0;
                                const LbmVec vel(mObjectSpeeds[OId]);
                                const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5;
                                USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
                                //errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
                                
                                for(size_t n=0; n<mMOIVertices.size(); n++) {
                                        POS2GRID_CHECK(mMOIVertices,n);
                                        // TODO - also reinit interface cells !?
                                        if(RFLAG(level, i,j,k, workSet)!=CFEmpty) { 
                                                // test prevent particle gen for inflow inter cells
                                                if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag
                                                continue; }
                                        monFluids++;

                                        // TODO add OPT3D treatment
                                        LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
                                        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;
                                        CellFlagType setFlag = CFInter;
                                        changeFlag(level, i,j,k, workSet, setFlag);
                                        mInitialMass += iniRho;
                                }
                                // second static init pass
                                if(staticInit) {
                                        CellFlagType set2Flag = CFMbndInflow|(OId<<24);
                                        for(size_t n=0; n<mMOIVertices.size(); n++) {
                                                POS2GRID_CHECK(mMOIVertices,n);
                                                if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
                                                if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
                                                        forceChangeFlag(level, i,j,k, workSet, set2Flag);
                                                } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
                                                        forceChangeFlag(level, i,j,k, workSet, 
                                                                        (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
                                                }
                                        }
                                } // second static inflow pass

                        } // inflow

                        else if(ntype&CFMbndOutflow){
                                const LbmFloat iniRho = 0.0;
                                for(size_t n=0; n<mMOIVertices.size(); n++) {
                                        POS2GRID_CHECK(mMOIVertices,n);
                                        // 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)) {
                                                        forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
                                                continue;
                                        }
                                        monFluids++;
                                        // interface cell - might be double empty? FIXME check?

                                        // remove fluid cells, shouldnt be here anyway
                                        LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
                                        LbmFloat fluidRho = RAC(tcel,0); FORDF1 { fluidRho += RAC(tcel,l); }
                                        mInitialMass -= fluidRho;
                                        RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
                                        RAC(tcel, dFlux) = FLUX_INIT;
                                        CellFlagType setFlag = CFInter;
                                        changeFlag(level, i,j,k, workSet, setFlag);

                                        // same as ifemptied for if below
                                        LbmPoint emptyp;
                                        emptyp.x = i; emptyp.y = j; emptyp.z = k;
                                        mListEmpty.push_back( emptyp );
                                        //calcCellsEmptied++;
                                } // points
                                // second static init pass
                                if(staticInit) {
                                        CellFlagType set2Flag = CFMbndOutflow|(OId<<24);
                                        for(size_t n=0; n<mMOIVertices.size(); n++) {
                                                POS2GRID_CHECK(mMOIVertices,n);
                                                if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
                                                if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
                                                        forceChangeFlag(level, i,j,k, workSet, set2Flag);
                                                } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
                                                        forceChangeFlag(level, i,j,k, workSet, 
                                                                        (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
                                                }
                                        }
                                } // second static outflow pass
                        } // outflow

                } // allbound check
        } // OId


        /* { // DEBUG check
                int workSet = mLevel[level].setCurr;
                FSGR_FORIJK1(level) {
                        if( (RFLAG(level,i,j,k, workSet) & CFBndMoving) ) {
                                if(QCELL(level, i,j,k, workSet, dFlux)!=targetTime) {
                                        errMsg("lastt"," old val!? at "<<PRINT_IJK<<" t="<<QCELL(level, i,j,k, workSet, dFlux)<<" target="<<targetTime);
                                }
                        }
                }
        } // DEBUG */
        
#undef POS2GRID_CHECK
        myTime_t monend = getTime();
        if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
        mLastSimTime = targetTime;
}


// 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); }

/*****************************************************************************/
/*****************************************************************************/
extern int globGeoInitDebug; //solver_interface
bool LbmFsgrSolver::initGeometryFlags() {
        int level = mMaxRefine;
        myTime_t geotimestart = getTime(); 
        ntlGeometryObject *pObj;
        ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
        debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
        // WARNING - copied to movobj init!

        /* init object velocities, this has always to be called for init */
        initGeoTree();
        if(mAllfluid) { 
                freeGeoTree();
                return true; }

        // make sure moving obstacles are inited correctly
        // preinit moving obj points
        int numobjs = (int)(mpGiObjects->size());
        for(int o=0; o<numobjs; o++) {
                ntlGeometryObject *obj = (*mpGiObjects)[o];
                //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
                if(
                                ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
                                (obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
                        if(!obj->getMeshAnimated()) { 
                                debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
                                obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
                        }
                }
        }

        // max speed init
        ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
        ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
        mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
        LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
        debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
        if(mpParam->getSimulationMaxSpeed() > allowMax) {
                // similar to adaptTimestep();
                LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
                LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
                debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
                mpParam->setDesiredTimestep( newdt );
                mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
                maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
                                      mSimulationTime, mpParam->getTimestep() )) ));
                debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
        }
        recalculateObjectSpeeds();
        // */

        // init obstacles for first time step (requires obj speeds)
        initMovingObstacles(true);

        /* set interface cells */
        ntlVec3Gfx pos,iniPos; // position of current cell
        LbmFloat rhomass = 0.0;
        CellFlagType ntype = CFInvalid;
        int savedNodes = 0;
        int OId = -1;
        gfxReal distance;

        // 2d display as rectangles
        if(LBMDIM==2) {
                dvec[2] = 0.0; 
                iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
                //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
        } else {
                iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
        }


        // first init boundary conditions
        // invalid cells are set to empty afterwards
        // TODO use floop macros!?
        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;
                                
                                GETPOS(i,j,k);
                                const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
                                if(inside) {
                                        pObj = (*mpGiObjects)[OId];
                                        switch( pObj->getGeoInitType() ){
                                        case FGI_MBNDINFLOW:  
                                          if(! pObj->getIsAnimated() ) {
                                                        rhomass = 1.0;
                                                        ntype = CFFluid | CFMbndInflow;
                                                } else {
                                                        ntype = CFInvalid;
                                                }
                                                break;
                                        case FGI_MBNDOUTFLOW: 
                                          if(! pObj->getIsAnimated() ) {
                                                        rhomass = 0.0;
                                                        ntype = CFEmpty|CFMbndOutflow; 
                                                } else {
                                                        ntype = CFInvalid;
                                                }
                                                break;
                                        case FGI_BNDNO: 
                                                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;
                                        default: // warn here?
                                                rhomass = BND_FILL;
                                                ntype = CFBnd|CFBndNoslip; break;
                                        }
                                }
                                if(ntype != CFInvalid) {
                                        // initDefaultCell
                                        if((ntype & CFMbndInflow) || (ntype & CFMbndOutflow) ) { }
                                        ntype |= (OId<<24); // NEWTEST2
                                        initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
                                }

                                // walk along x until hit for following inits
                                if(distance<=-1.0) { distance = 100.0; } // FIXME dangerous
                                if(distance>=0.0) {
                                        gfxReal dcnt=dvec[0];
                                        while(( dcnt< distance-dvec[0] )&&(i+1<mLevel[level].lSizex-1)) {
                                                dcnt += dvec[0]; i++;
                                                savedNodes++;
                                                if(ntype != CFInvalid) {
                                                        // rho,mass,OId are still inited from above
                                                        initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
                                                }
                                        }
                                } 
                                // */

                        } 
                } 
        } // zmax
        // */

        // now init fluid layer
        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++) {
                                if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
                                ntype = CFInvalid;
                                int inits = 0;
                                GETPOS(i,j,k);
                                const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
                                if(inside) {
                                        ntype = CFFluid;
                                }
                                if(ntype != CFInvalid) {
                                        // initDefaultCell
                                        rhomass = 1.0;
                                        initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
                                        inits++;
                                }

                                // walk along x until hit for following inits
                                if(distance<=-1.0) { distance = 100.0; }
                                if(distance>=0.0) {
                                        gfxReal dcnt=dvec[0];
                                        while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
                                                dcnt += dvec[0]; i++;
                                                savedNodes++;
                                                if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
                                                if(ntype != CFInvalid) {
                                                        // rhomass are still inited from above
                                                        initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
                                                        inits++;
                                                }
                                        }
                                } // distance>0
                                
                        } 
                } 
        } // zmax

        // reset invalid to empty again
        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++) {
                                if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) {
                                        RFLAG(level, i,j,k, mLevel[level].setOther) = 
                                        RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty;
                                }
                        }
                }
        }

        freeGeoTree();
        myTime_t geotimeend = getTime(); 
        debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 ); 
        //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
        return true;
}

#undef GETPOS

/*****************************************************************************/
/* init part for all freesurface testcases */
void LbmFsgrSolver::initFreeSurfaces() {
        double interfaceFill = 0.45;   // filling level of interface cells
        //interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!

        // set interface cells 
        FSGR_FORIJK1(mMaxRefine) {
                if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
                        FORDF1 {
                                int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
                                if( ( 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) );
                                        // unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
                                        initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
                                }
                        }
                }
        }

        // remove invalid interface cells 
        FSGR_FORIJK1(mMaxRefine) {
                // remove invalid interface cells 
                if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
                        int delit = 0;
                        int NBs = 0; // neighbor flags or'ed 
                        int noEmptyNB = 1;

                        FORDF1 {
                                if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
                                        noEmptyNB = 0;
                                }
                                NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
                        }
                        // remove cells with no fluid or interface neighbors
                        if((NBs & CFFluid)==0) { delit = 1; }
                        if((NBs & CFInter)==0) { delit = 1; }
                        // remove cells with no empty neighbors
                        if(noEmptyNB) { delit = 2; }

                        // now we can remove the cell 
                        if(delit==1) {
                                initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
                        }
                        if(delit==2) {
                                //initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
                                LbmFloat arho=0., aux=0., auy=0., auz=0.;
                                interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
                                initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
                        }
                } // interface 
        } // */

        // another brute force init, make sure the fill values are right...
        // and make sure both sets are equal
        for(int lev=0; lev<=mMaxRefine; lev++) {
        FSGR_FORIJK_BOUNDS(lev) {
                if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { 
                        QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
                        continue;
                }
                if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { 
                        QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
                        continue;
                }
        } }

        // ----------------------------------------------------------------------
        // smoother surface...
        if(mInitSurfaceSmoothing>0) {
                debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing timestep ",10);
#if COMPRESSGRIDS==1
                //errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
#endif // COMPRESSGRIDS==0
        }
        for(int s=0; s<mInitSurfaceSmoothing; s++) {
                //SGR_FORIJK1(mMaxRefine) {

                int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
                int lev = 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=1;j<mLevel[lev].lSizey-1;++j) {
                for(int i=1;i<mLevel[lev].lSizex-1;++i) {
                        if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
                                LbmFloat mass = 0.0;
                                //LbmFloat nbdiv;
                                //FORDF0 {
                                for(int l=0;(l<cDirNum); l++) { 
                                        int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
                                        if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
                                                mass += 1.0;
                                        }
                                        if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFInter ){
                                                mass += QCELL(lev, ni,nj,nk, mLevel[lev].setCurr, dMass);
                                        }
                                        //nbdiv+=1.0;
                                }

                                //errMsg(" I ", PRINT_IJK<<" m"<<mass );
                                QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
                                QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
                        }
                }}}

                mLevel[lev].setOther = mLevel[lev].setCurr;
                mLevel[lev].setCurr ^= 1;
        }
        // copy back...?
}

/*****************************************************************************/
/* init part for all freesurface testcases */
void LbmFsgrSolver::initStandingFluidGradient() {
        // ----------------------------------------------------------------------
        // standing fluid preinit
        const int debugStandingPreinit = 0;
        int haveStandingFluid = 0;

#define STANDFLAGCHECK(iindex) \
                                if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
                                                ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
                                        if((iindex)>1) { haveStandingFluid=(iindex); } \
                                        j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
                                        continue; \
                                } 
        int gravIndex[3] = {0,0,0};
        int gravDir[3] = {1,1,1};
        int maxGravComp = 1; // by default y
        int gravComp1 = 0; // by default y
        int gravComp2 = 2; // by default y
        if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
        if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }

        int gravIMin[3] = { 0 , 0 , 0 };
        int gravIMax[3] = {
                mLevel[mMaxRefine].lSizex + 0,
                mLevel[mMaxRefine].lSizey + 0,
                mLevel[mMaxRefine].lSizez + 0 };
        if(LBMDIM==2) gravIMax[2] = 1;

        //int gravDir = 1;
        if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
                // swap directions
                int i=maxGravComp;
                int tmp = gravIMin[i];
                gravIMin[i] = gravIMax[i] - 1;
                gravIMax[i] = tmp - 1;
                gravDir[i] = -1;
        }
#define PRINTGDIRS \
        errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
        errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
        errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); 
        // _PRINTGDIRS;

        bool gravAbort = false;
#define GRAVLOOP \
        gravAbort=false; \
        for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
                for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
                        for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0]) 

        GRAVLOOP {
                int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
                if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || 
                    ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) || 
                                ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) {  
                        int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
                        if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
                        if(fluidHeight>1) {
                                haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; 
                                gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
                        }
                        gravAbort = true; continue; 
                } 
        } // GRAVLOOP
        // _PRINTGDIRS;

        LbmFloat fluidHeight;
        fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
#if LBM_INCLUDE_TESTSOLVERS==1
        if(mUseTestdata) mpTest->mFluidHeight = (int)fluidHeight;
#endif // ELBEEM_PLUGIN!=1
        if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "fheight="<<fluidHeight<<" min="<<PRINT_VEC(gravIMin[0],gravIMin[1],       gravIMin[2])<<" max="<<PRINT_VEC(gravIMax[0], gravIMax[1],gravIMax[2])<<
                        " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
                                
        if(mDisableStandingFluidInit) {
                debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
                haveStandingFluid=0;
        }

        // copy flags and init , as no flags will be changed during grav init
        // also important for corasening later on
        const int lev = mMaxRefine;
        CellFlagType nbflag[LBM_DFNUM], nbored; 
        for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
                for(int j=0;j<mLevel[lev].lSizey-0;++j) {
                        for(int i=0;i<mLevel[lev].lSizex-0;++i) {
                                if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
                                        nbored = 0;
                                        FORDF1 {
                                                nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
                                                nbored |= nbflag[l];
                                        } 
                                        if(nbored&CFBnd) {
                                                RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
                                        } else {
                                                RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
                                        }
                                }
                                RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
        } } }

        if(haveStandingFluid) {
                int rhoworkSet = mLevel[lev].setCurr;
                myTime_t timestart = getTime(); // FIXME use user time here?

                GRAVLOOP {
                        int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
                        //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
                        if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
                                        ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ 
                                //gravAbort = true; 
                                continue;
                        }

                        LbmFloat rho = 1.0;
                        // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
                        rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * 
                                (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
                        if(debugStandingPreinit) 
                                if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { 
                                        errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); 
                                }

                        if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
                                        (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
                                FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
                                QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
                        }

                } // 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);
                
                int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
                preinitSteps = (haveStandingFluid>>2); // not much use...?
                //preinitSteps = 0;
                debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
                for(int s=0; s<preinitSteps; s++) {
                        // in solver main cpp
                        standingFluidPreinit();
                }

                myTime_t timeend = getTime();
                debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
#undef  NBFLAG
        }
}



bool LbmFsgrSolver::checkSymmetry(string idstring)
{
        bool erro = false;
        bool symm = true;
        int msgs = 0;
        const int maxMsgs = 10;
        const bool markCells = false;

        //for(int lev=0; lev<=mMaxRefine; lev++) {
        { int lev = mMaxRefine;

        // no point if not symm.
        if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
                // ok
        } else {
                return false;
        }

        for(int s=0; s<2; s++) {
        FSGR_FORIJK1(lev) {
                if(i<(mLevel[lev].lSizex/2)) {
                        int inb = (mLevel[lev].lSizey-1-i); 

                        if(lev==mMaxRefine) inb -= 1;           // FSGR_SYMM_T

                        if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
                                if(LBMDIM==2) {
                                        if(msgs<maxMsgs) { msgs++;
                                                errMsg("EFLAG", PRINT_IJK<<"s"<<s<<" flag "<<RFLAG(lev, i,j,k,s)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" flag "<<RFLAG(lev, inb,j,k,s) );
                                        }
                                }
                                if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
                                symm = false;
                        }
                        if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
                                if(LBMDIM==2) {
                                        if(msgs<maxMsgs) { msgs++;
                                                //debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
                                                errMsg("EMASS", PRINT_IJK<<"s"<<s<<" mass "<<QCELL(lev, i,j,k,s, dMass)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" mass "<<QCELL(lev, inb,j,k,s, dMass) );
                                        }
                                }
                                if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
                                symm = false;
                        }

                        LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
                        FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
                        LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
                        FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
                        if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
                                if(LBMDIM==2) {
                                        if(msgs<maxMsgs) { msgs++;
                                                //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
                                                errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
                                        }
                                }
                                if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
                                symm = false;
                        }
                }
        } }
        } // lev
        LbmFloat maxdiv =0.0;
        if(erro) {
                errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
                //if(LBMDIM==2) mPanic = true; 
                //return false;
        } else {
                errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
        }
        // all ok...
        return symm;
}// */


void 
LbmFsgrSolver::interpolateCellValues(
                int level,int ei,int ej,int ek,int workSet, 
                LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz) 
{
        LbmFloat avgrho = 0.0;
        LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
        LbmFloat cellcnt = 0.0;

        for(int nbl=1; nbl< cDfNum ; ++nbl) {
                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< cDfNum ; ++rl) { 
                                LbmFloat nbdf =  QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
                                avgux  += (dfDvecX[rl]*nbdf); 
                                avguy  += (dfDvecY[rl]*nbdf);  
                                avguz  += (dfDvecZ[rl]*nbdf);  
                                avgrho += nbdf;
                        }
                }
        }

        if(cellcnt<=0.0) {
                // no nbs? just use eq.
                avgrho = 1.0;
                avgux = avguy = avguz = 0.0;
                //TTT mNumProblems++;
#if ELBEEM_PLUGIN!=1
                //mPanic=1; 
                // this can happen for worst case moving obj scenarios...
                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;
        }

        retrho = avgrho;
        retux = avgux;
        retuy = avguy;
        retuz = avguz;
} // interpolateCellValues


/******************************************************************************
 * instantiation
 *****************************************************************************/

LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }