intern/elbeem/intern/solver_init.cpp

   1 /******************************************************************************
   2  *
   3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
   4  * Copyright 2003-2006 Nils Thuerey
   5  *
   6  * Standard LBM Factory implementation
   7  *
   8  *****************************************************************************/
   9
  10
  11 #include "solver_class.h"
  12 #include "solver_relax.h"
  13 // for geo init FGI_ defines
  14 #include "elbeem.h"
  15
  16 // helper for 2d init
  17 #define SWAPYZ(vec) { \
  18                 const LbmFloat tmp = (vec)[2]; \
  19                 (vec)[2] = (vec)[1]; (vec)[1] = tmp; }
  20
  21
  22 /*****************************************************************************/
  23 //! common variables
  24
  25 /*****************************************************************************/
  26 /*! 3D implementation D3Q19 */
  27 #if LBMDIM==3
  28
  29         //! how many dimensions?
  30         const int LbmFsgrSolver::cDimension = 3;
  31
  32         // Wi factors for collide step
  33         const LbmFloat LbmFsgrSolver::cCollenZero    = (1.0/3.0);
  34         const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/18.0);
  35         const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
  36
  37         //! threshold value for filled/emptied cells
  38         const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
  39         const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
  40         const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
  41         const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;
  42
  43         //! size of a single set of distribution functions
  44         const int    LbmFsgrSolver::cDfNum      = 19;
  45         //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
  46         const int    LbmFsgrSolver::cDirNum     = 27;
  47
  48         //const string LbmFsgrSolver::dfString[ cDfNum ] = {
  49         const char* LbmFsgrSolver::dfString[ cDfNum ] = {
  50                 " C", " N"," S"," E"," W"," T"," B",
  51                 "NE","NW","SE","SW",
  52                 "NT","NB","ST","SB",
  53                 "ET","EB","WT","WB"
  54         };
  55
  56         const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
  57                 cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB,
  58                 cDirNE, cDirNW, cDirSE, cDirSW,
  59                 cDirNT, cDirNB, cDirST, cDirSB,
  60                 cDirET, cDirEB, cDirWT, cDirWB
  61         };
  62
  63         const int LbmFsgrSolver::dfInv[ cDfNum ] = {
  64                 cDirC,  cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
  65                 cDirSW, cDirSE, cDirNW, cDirNE,
  66                 cDirSB, cDirST, cDirNB, cDirNT,
  67                 cDirWB, cDirWT, cDirEB, cDirET
  68         };
  69
  70         const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
  71                 0,  0, 0, 0, 0, 0, 0,
  72                 cDirSE, cDirSW, cDirNE, cDirNW,
  73                 0, 0, 0, 0,
  74                 cDirEB, cDirET, cDirWB, cDirWT
  75         };
  76
  77         const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
  78                 0,  0, 0, 0, 0, 0, 0,
  79                 cDirNW, cDirNE, cDirSW, cDirSE,
  80                 cDirNB, cDirNT, cDirSB, cDirST,
  81                 0, 0, 0, 0
  82         };
  83
  84         const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
  85                 0,  0, 0, 0, 0, 0, 0,
  86                 0, 0, 0, 0,
  87                 cDirST, cDirSB, cDirNT, cDirNB,
  88                 cDirWT, cDirWB, cDirET, cDirEB
  89         };
  90
  91         // Vector Order 3D:
  92         //  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
  93         //  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
  94         //  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
  95         //  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
  96
  97         const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
  98                 0, 0,0, 1,-1, 0,0,
  99                 1,-1,1,-1,
 100                 0,0,0,0,
 101                 1,1,-1,-1,
 102                  1,-1, 1,-1,
 103                  1,-1, 1,-1,
 104         };
 105         const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
 106                 0, 1,-1, 0,0,0,0,
 107                 1,1,-1,-1,
 108                 1,1,-1,-1,
 109                 0,0,0,0,
 110                  1, 1,-1,-1,
 111                  1, 1,-1,-1
 112         };
 113         const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
 114                 0, 0,0,0,0,1,-1,
 115                 0,0,0,0,
 116                 1,-1,1,-1,
 117                 1,-1,1,-1,
 118                  1, 1, 1, 1,
 119                 -1,-1,-1,-1
 120         };
 121
 122         const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
 123                 0, 0,0, 1,-1, 0,0,
 124                 1,-1,1,-1,
 125                 0,0,0,0,
 126                 1,1,-1,-1,
 127                  1,-1, 1,-1,
 128                  1,-1, 1,-1
 129         };
 130         const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
 131                 0, 1,-1, 0,0,0,0,
 132                 1,1,-1,-1,
 133                 1,1,-1,-1,
 134                 0,0,0,0,
 135                  1, 1,-1,-1,
 136                  1, 1,-1,-1
 137         };
 138         const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
 139                 0, 0,0,0,0,1,-1,
 140                 0,0,0,0,
 141                 1,-1,1,-1,
 142                 1,-1,1,-1,
 143                  1, 1, 1, 1,
 144                 -1,-1,-1,-1
 145         };
 146
 147         /* principal directions */
 148         const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
 149                 1,-1, 0,0, 0,0
 150         };
 151         const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
 152                 0,0, 1,-1, 0,0
 153         };
 154         const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
 155                 0,0, 0,0, 1,-1
 156         };
 157
 158         /*! arrays for les model coefficients, inited in lbmsolver constructor */
 159         LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
 160         LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
 161
 162
 163         const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
 164                 cCollenZero,
 165                 cCollenOne, cCollenOne, cCollenOne,
 166                 cCollenOne, cCollenOne, cCollenOne,
 167                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
 168                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
 169                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
 170         };
 171
 172         /* precalculated equilibrium dfs, inited in lbmsolver constructor */
 173         LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
 174
 175 #else // end LBMDIM==3 , LBMDIM==2
 176
 177 /*****************************************************************************/
 178 /*! 2D implementation D2Q9 */
 179
 180                 //! how many dimensions?
 181                 const int LbmFsgrSolver::cDimension = 2;
 182
 183                 //! Wi factors for collide step
 184                 const LbmFloat LbmFsgrSolver::cCollenZero    = (4.0/9.0);
 185                 const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/9.0);
 186                 const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
 187
 188                 //! threshold value for filled/emptied cells
 189                 const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
 190                 const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
 191                 const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
 192                 const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;
 193
 194                 //! size of a single set of distribution functions
 195                 const int LbmFsgrSolver::cDfNum  = 9;
 196                 const int LbmFsgrSolver::cDirNum = 9;
 197
 198         //const string LbmFsgrSolver::dfString[ cDfNum ] = {
 199         const char* LbmFsgrSolver::dfString[ cDfNum ] = {
 200                 " C",
 201                 " N",   " S", " E", " W",
 202                 "NE", "NW", "SE","SW"
 203         };
 204
 205         const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
 206                 cDirC,
 207                 cDirN,  cDirS,  cDirE,  cDirW,
 208                 cDirNE, cDirNW, cDirSE, cDirSW
 209         };
 210
 211         const int LbmFsgrSolver::dfInv[ cDfNum ] = {
 212                 cDirC,
 213                 cDirS,  cDirN,  cDirW,  cDirE,
 214                 cDirSW, cDirSE, cDirNW, cDirNE
 215         };
 216
 217         const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
 218                 0,
 219                 0,  0,  0,  0,
 220                 cDirSE, cDirSW, cDirNE, cDirNW
 221         };
 222
 223         const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
 224                 0,
 225                 0,  0,  0,  0,
 226                 cDirNW, cDirNE, cDirSW, cDirSE
 227         };
 228
 229         const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
 230                 0,  0, 0, 0, 0,
 231                 0, 0, 0, 0
 232         };
 233
 234         // Vector Order 2D:
 235         // 0  1 2  3  4  5  6 7  8
 236         // 0, 0,0, 1,-1, 1,-1,1,-1
 237         // 0, 1,-1, 0,0, 1,1,-1,-1
 238
 239         const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
 240                 0,
 241                 0,0, 1,-1,
 242                 1,-1,1,-1
 243         };
 244         const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
 245                 0,
 246                 1,-1, 0,0,
 247                 1,1,-1,-1
 248         };
 249         const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
 250                 0, 0,0,0,0, 0,0,0,0
 251         };
 252
 253         const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
 254                 0,
 255                 0,0, 1,-1,
 256                 1,-1,1,-1
 257         };
 258         const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
 259                 0,
 260                 1,-1, 0,0,
 261                 1,1,-1,-1
 262         };
 263         const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
 264                 0, 0,0,0,0, 0,0,0,0
 265         };
 266
 267         const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
 268                 1,-1, 0,0
 269         };
 270         const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
 271                 0,0, 1,-1
 272         };
 273         const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
 274                 0,0, 0,0
 275         };
 276
 277
 278         /*! arrays for les model coefficients, inited in lbmsolver constructor */
 279         LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
 280         LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
 281
 282
 283         const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
 284                 cCollenZero,
 285                 cCollenOne, cCollenOne, cCollenOne, cCollenOne,
 286                 cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
 287         };
 288
 289         /* precalculated equilibrium dfs, inited in lbmsolver constructor */
 290         LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
 291
 292 // D2Q9 end
 293 #endif  // LBMDIM==2
 294
 295
 296 // required globals
 297 extern bool glob_mpactive;
 298 extern int glob_mpnum, glob_mpindex;
 299
 300
 301 /******************************************************************************
 302  * Lbm Constructor
 303  *****************************************************************************/
 304 LbmFsgrSolver::LbmFsgrSolver() :
 305         //D(),
 306         mCurrentMass(0.0), mCurrentVolume(0.0),
 307         mNumProblems(0),
 308         mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
 309         mpPreviewSurface(NULL),
 310         mTimeAdap(true), mForceTimeStepReduce(false),
 311         mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
 312         mInitSurfaceSmoothing(0), mFsSurfGenSetting(0),
 313         mTimestepReduceLock(0),
 314         mTimeSwitchCounts(0), mTimeMaxvelStepCnt(0),
 315         mSimulationTime(0.0), mLastSimTime(0.0),
 316         mMinTimestep(0.0), mMaxTimestep(0.0),
 317         mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
 318         mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
 319         mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
 320         mIsoWeightMethod(1),
 321         mMaxRefine(1),
 322         mDfScaleUp(-1.0), mDfScaleDown(-1.0),
 323         mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03
 324         mDebugOmegaRet(0.0),
 325         mLastOmega(1e10), mLastGravity(1e10),
 326         mNumInvIfTotal(0), mNumFsgrChanges(0),
 327         mDisableStandingFluidInit(0),
 328         mInit2dYZ(false),
 329         mForceTadapRefine(-1), mCutoff(-1)
 330 {
 331 #if LBM_INCLUDE_CONTROL==1
 332         mpControl = new LbmControlData();
 333 #endif
 334
 335 #if LBM_INCLUDE_TESTSOLVERS==1
 336         mpTest = new LbmTestdata();
 337         mMpNum = mMpIndex = 0;
 338         mOrgSizeX = 0;
 339         mOrgStartX = 0.;
 340         mOrgEndX = 0.;
 341 #endif // LBM_INCLUDE_TESTSOLVERS!=1
 342         mpIso = new IsoSurface( mIsoValue );
 343
 344   // init equilibrium dist. func
 345   LbmFloat rho=1.0;
 346   FORDF0 {
 347                 dfEquil[l] = this->getCollideEq( l,rho,  0.0, 0.0, 0.0);
 348   }
 349         dfEquil[dMass] = 1.;
 350         dfEquil[dFfrac] = 1.;
 351         dfEquil[dFlux] = FLUX_INIT;
 352
 353         // init LES
 354         int odm = 0;
 355         for(int m=0; m<LBMDIM; m++) {
 356                 for(int l=0; l<cDfNum; l++) {
 357                         this->lesCoeffDiag[m][l] =
 358                         this->lesCoeffOffdiag[m][l] = 0.0;
 359                 }
 360         }
 361         for(int m=0; m<LBMDIM; m++) {
 362                 for(int n=0; n<LBMDIM; n++) {
 363                         for(int l=1; l<cDfNum; l++) {
 364                                 LbmFloat em;
 365                                 switch(m) {
 366                                         case 0: em = dfDvecX[l]; break;
 367                                         case 1: em = dfDvecY[l]; break;
 368                                         case 2: em = dfDvecZ[l]; break;
 369                                         default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
 370                                 }
 371                                 LbmFloat en;
 372                                 switch(n) {
 373                                         case 0: en = dfDvecX[l]; break;
 374                                         case 1: en = dfDvecY[l]; break;
 375                                         case 2: en = dfDvecZ[l]; break;
 376                                         default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
 377                                 }
 378                                 const LbmFloat coeff = em*en;
 379                                 if(m==n) {
 380                                         this->lesCoeffDiag[m][l] = coeff;
 381                                 } else {
 382                                         if(m>n) {
 383                                                 this->lesCoeffOffdiag[odm][l] = coeff;
 384                                         }
 385                                 }
 386                         }
 387
 388                         if(m==n) {
 389                         } else {
 390                                 if(m>n) odm++;
 391                         }
 392                 }
 393         }
 394
 395         mDvecNrm[0] = LbmVec(0.0);
 396   FORDF1 {
 397                 mDvecNrm[l] = getNormalized(
 398                         LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] )
 399                         ) * -1.0;
 400         }
 401
 402         // calculate gauss weights for restriction
 403         //LbmFloat mGaussw[27];
 404         LbmFloat totGaussw = 0.0;
 405         const LbmFloat alpha = 1.0;
 406         const LbmFloat gw = sqrt(2.0*LBMDIM);
 407 #if ELBEEM_PLUGIN!=1
 408         errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
 409 #endif
 410         for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
 411         //for(int n=0;(n<cDirNum); n++) {
 412         for(int n=0;(n<cDfNum); n++) {
 413                 const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
 414                 LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
 415                 mGaussw[n] = w;
 416                 totGaussw += w;
 417         }
 418         for(int n=0;(n<cDirNum); n++) {
 419                 mGaussw[n] = mGaussw[n]/totGaussw;
 420         }
 421
 422 }
 423
 424 /*****************************************************************************/
 425 /* Destructor */
 426 /*****************************************************************************/
 427 LbmFsgrSolver::~LbmFsgrSolver()
 428 {
 429   if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; }
 430 #if COMPRESSGRIDS==1
 431         delete [] mLevel[mMaxRefine].mprsCells[1];
 432         mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
 433 #endif // COMPRESSGRIDS==1
 434
 435         for(int i=0; i<=mMaxRefine; i++) {
 436                 for(int s=0; s<2; s++) {
 437                         if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
 438                         if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
 439                 }
 440         }
 441         delete mpIso;
 442         if(mpPreviewSurface) delete mpPreviewSurface;
 443         // cleanup done during scene deletion...
 444
 445         // always output performance estimate
 446         debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
 447   if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
 448 }
 449
 450
 451
 452
 453 /******************************************************************************
 454  * initilize variables fom attribute list
 455  *****************************************************************************/
 456 void LbmFsgrSolver::parseAttrList()
 457 {
 458         LbmSolverInterface::parseStdAttrList();
 459
 460         string matIso("default");
 461         matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
 462         mpIso->setMaterialName( matIso );
 463         mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
 464         mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
 465         mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false );
 466         mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );
 467
 468         mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
 469         mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
 470         mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
 471         mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
 472         mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
 473
 474         // refinement
 475         mMaxRefine = mRefinementDesired;
 476         mMaxRefine  = mpSifAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
 477         if(mMaxRefine<0) mMaxRefine=0;
 478         if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
 479         mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
 480         mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
 481         mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
 482
 483         // demo mode settings
 484         mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
 485         // FIXME check needed?
 486         mFVArea   = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
 487
 488         // debugging - skip some time...
 489         double starttimeskip = 0.;
 490         starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
 491         mSimulationTime += starttimeskip;
 492         if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
 493
 494 #if LBM_INCLUDE_CONTROL==1
 495         mpControl->parseControldataAttrList(mpSifAttrs);
 496 #endif
 497
 498 #if LBM_INCLUDE_TESTSOLVERS==1
 499         mUseTestdata = 0;
 500         mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
 501         mpTest->parseTestdataAttrList(mpSifAttrs);
 502 #ifdef ELBEEM_PLUGIN
 503         mUseTestdata=1; // DEBUG
 504 #endif // ELBEEM_PLUGIN
 505
 506         mMpNum = mpSifAttrs->readInt("mpnum",  mMpNum ,"LbmFsgrSolver", "mMpNum", false);
 507         mMpIndex = mpSifAttrs->readInt("mpindex",  mMpIndex ,"LbmFsgrSolver", "mMpIndex", false);
 508         if(glob_mpactive) {
 509                 // used instead...
 510                 mMpNum = glob_mpnum;
 511                 mMpIndex = glob_mpindex;
 512         } else {
 513                 glob_mpnum = mMpNum;
 514                 glob_mpindex = 0;
 515         }
 516         errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<<glob_mpactive<<", "<<glob_mpindex<<"/"<<glob_mpnum);
 517         if(mMpNum>0) {
 518                 mUseTestdata=1; // needed in this case...
 519         }
 520
 521         errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
 522 #else // LBM_INCLUDE_TESTSOLVERS!=1
 523         // not testsolvers, off by default
 524         mUseTestdata = 0;
 525         if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
 526 #endif // LBM_INCLUDE_TESTSOLVERS!=1
 527
 528         mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
 529         // deprecated!
 530         float mInitialCsmagoCoarse = 0.0;
 531         mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
 532 #if USE_LES==1
 533 #else // USE_LES==1
 534         debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
 535         mInitialCsmago = 0.0;
 536 #endif // USE_LES==1
 537 }
 538
 539
 540 /******************************************************************************
 541  * Initialize omegas and forces on all levels (for init/timestep change)
 542  *****************************************************************************/
 543 void LbmFsgrSolver::initLevelOmegas()
 544 {
 545         // no explicit settings
 546         mOmega = mpParam->calculateOmega(mSimulationTime);
 547         mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) );
 548         mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused
 549         if(mInit2dYZ) { SWAPYZ(mGravity); }
 550
 551         // check if last init was ok
 552         LbmFloat gravDelta = norm(mGravity-mLastGravity);
 553         //errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<mOmega<<" - lom:"<<mLastOmega<<" gv:"<<mGravity<<" - "<<mLastGravity<<" , "<<gravDelta  );
 554         if((mOmega == mLastOmega) && (gravDelta<=0.0)) return;
 555
 556         if(mInitialCsmago<=0.0) {
 557                 if(OPT3D==1) {
 558                         errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR);
 559                         return;
 560                 }
 561         }
 562
 563         LbmFloat  fineCsmago  = mInitialCsmago;
 564         LbmFloat coarseCsmago = mInitialCsmago;
 565         LbmFloat maxFineCsmago1    = 0.026;
 566         LbmFloat maxCoarseCsmago1  = 0.029; // try stabilizing
 567         LbmFloat maxFineCsmago2    = 0.028;
 568         LbmFloat maxCoarseCsmago2  = 0.032; // try stabilizing some more
 569         if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) {
 570                 fineCsmago = maxFineCsmago1;
 571                 coarseCsmago = maxCoarseCsmago1;
 572         }
 573         if((mMaxRefine>1)&&(mInitialCsmago<maxFineCsmago2)) {
 574                 fineCsmago = maxFineCsmago2;
 575                 coarseCsmago = maxCoarseCsmago2;
 576         }
 577
 578
 579         // use Tau instead of Omega for calculations
 580         { // init base level
 581                 int i = mMaxRefine;
 582                 mLevel[i].omega    = mOmega;
 583                 mLevel[i].timestep = mpParam->getTimestep();
 584                 mLevel[i].lcsmago = fineCsmago; //CSMAGO_INITIAL;
 585                 mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
 586                 mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
 587         }
 588
 589         // init all sub levels
 590         for(int i=mMaxRefine-1; i>=0; i--) {
 591                 //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
 592                 double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
 593                 nomega                = 1.0/nomega;
 594                 mLevel[i].omega       = (LbmFloat)nomega;
 595                 mLevel[i].timestep    = 2.0 * mLevel[i+1].timestep;
 596                 mLevel[i].lcsmago = coarseCsmago;
 597                 mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
 598                 mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
 599         }
 600
 601         // for lbgk
 602         mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
 603         for(int i=mMaxRefine-1; i>=0; i--) {
 604                 // should be the same on all levels...
 605                 // for lbgk
 606                 mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
 607         }
 608
 609         mLastOmega = mOmega;
 610         mLastGravity = mGravity;
 611         // debug? invalidate old values...
 612         mGravity = -100.0;
 613         mOmega = -100.0;
 614
 615         for(int i=0; i<=mMaxRefine; i++) {
 616                 if(!mSilent) {
 617                         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
 618                                         <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
 619                                         <<" cmsagp:"<<mLevel[i].lcsmago<<", "
 620                                         << " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
 621                 } else {
 622                         if(!mInitDone) {
 623                                 debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
 624                                                 <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
 625                         }
 626                 }
 627         }
 628         if(mMaxRefine>0) {
 629                 mDfScaleUp   = (mLevel[0  ].timestep/mLevel[0+1].timestep)* (1.0/mLevel[0  ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
 630                 mDfScaleDown = (mLevel[0+1].timestep/mLevel[0  ].timestep)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0  ].omega-1.0); // yu
 631         }
 632 }
 633
 634
 635 /******************************************************************************
 636  * Init Solver (values should be read from config file)
 637  *****************************************************************************/
 638
 639 /*! finish the init with config file values (allocate arrays...) */
 640 bool LbmFsgrSolver::initializeSolverMemory()
 641 {
 642   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<mInitDone<<" "<<(void*)this,1);
 643
 644         // init cppf stage
 645         if(mCppfStage>0) {
 646                 mSizex *= mCppfStage;
 647                 mSizey *= mCppfStage;
 648                 mSizez *= mCppfStage;
 649         }
 650         if(mFsSurfGenSetting==-1) {
 651                 // all on
 652                 mFsSurfGenSetting =
 653                          fssgNormal   | fssgNoNorth  | fssgNoSouth  | fssgNoEast   |
 654                          fssgNoWest   | fssgNoTop    | fssgNoBottom | fssgNoObs   ;
 655         }
 656
 657         // size inits to force cubic cells and mult4 level dimensions
 658         // and make sure we dont allocate too much...
 659         bool memOk = false;
 660         int orgSx = mSizex;
 661         int orgSy = mSizey;
 662         int orgSz = mSizez;
 663         double sizeReduction = 1.0;
 664         double memEstFromFunc = -1.0;
 665         double memEstFine = -1.0;
 666         string memreqStr("");
 667         bool firstMInit = true;
 668         int minitTries=0;
 669         while(!memOk) {
 670                 minitTries++;
 671                 initGridSizes( mSizex, mSizey, mSizez,
 672                                 mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL);
 673
 674                 // MPT
 675 #if LBM_INCLUDE_TESTSOLVERS==1
 676                 if(firstMInit) {
 677                         mrSetup();
 678                 }
 679 #endif // LBM_INCLUDE_TESTSOLVERS==1
 680                 firstMInit=false;
 681
 682                 calculateMemreqEstimate( mSizex, mSizey, mSizez,
 683                                 mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr );
 684
 685                 double memLimit;
 686                 string memLimStr("-");
 687                 if(sizeof(void*)==4) {
 688                         // 32bit system, limit to 2GB
 689                         memLimit = 2.0* 1024.0*1024.0*1024.0;
 690                         memLimStr = string("2GB");
 691                 } else {
 692                         // 64bit, just take 16GB as limit for now...
 693                         memLimit = 16.0* 1024.0*1024.0*1024.0;
 694                         memLimStr = string("16GB");
 695                 }
 696
 697                 // restrict max. chunk of 1 mem block to 1GB for windos
 698                 bool memBlockAllocProblem = false;
 699                 double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
 700                 //std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
 701 #ifdef WIN32
 702                 double maxWinMemChunk = 1100.*1024.*1024.;
 703                 if(memEstFine> maxWinMemChunk) {
 704                         memBlockAllocProblem = true;
 705                 }
 706 #endif // WIN32
 707 #ifdef __APPLE__
 708                 double maxMacMemChunk = 1200.*1024.*1024.;
 709                 if(memEstFine> maxMacMemChunk) {
 710                         memBlockAllocProblem = true;
 711                 }
 712 #endif // Mac
 713                 if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
 714                         // max memory chunk for 32bit systems 2gig
 715                         memBlockAllocProblem = true;
 716                 }
 717
 718                 if(memEstFromFunc>memLimit || memBlockAllocProblem) {
 719                         sizeReduction *= 0.9;
 720                         mSizex = (int)(orgSx * sizeReduction);
 721                         mSizey = (int)(orgSy * sizeReduction);
 722                         mSizez = (int)(orgSz * sizeReduction);
 723                         debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
 724                                         //memEstFromFunc<<"/"<<memLimit<<", "<<
 725                                         //memEstFine<<"/"<<maxWinMemChunk<<", "<<
 726                                         memreqStr<<"/"<<memLimStr<<", "<<
 727                                         "retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
 728                                         , 3 );
 729                 } else {
 730                         memOk = true;
 731                 }
 732         }
 733
 734         mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
 735   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
 736           ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
 737           ", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
 738           ", est. Mem.Req.: "<<memreqStr        ,2);
 739         mpParam->setSize(mSizex, mSizey, mSizez);
 740         if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }
 741
 742   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
 743                 <<"LBM_EPSILON="<<LBM_EPSILON  <<" "
 744                 <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG <<" "
 745                 <<"OPT3D="<<OPT3D <<" "
 746                 <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
 747                 <<"MASS_INVALID="<<MASS_INVALID <<" "
 748                 <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
 749                 <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
 750                 <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
 751                 <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" "
 752                 <<"USE_LES="<<USE_LES <<" "
 753                 ,10);
 754
 755         // perform 2D corrections...
 756         if(LBMDIM == 2) mSizez = 1;
 757
 758         mpParam->setSimulationMaxSpeed(0.0);
 759         if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
 760         mpParam->setTadapLevels( mMaxRefine+1 );
 761
 762         if(mForceTadapRefine>mMaxRefine) {
 763                 mpParam->setTadapLevels( mForceTadapRefine+1 );
 764                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
 765         }
 766
 767         if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
 768                 errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
 769                 return false;
 770         }
 771
 772
 773         // init vectors
 774         for(int i=0; i<=mMaxRefine; i++) {
 775                 mLevel[i].id = i;
 776                 mLevel[i].nodeSize = 0.0;
 777                 mLevel[i].simCellSize = 0.0;
 778                 mLevel[i].omega = 0.0;
 779                 mLevel[i].time = 0.0;
 780                 mLevel[i].timestep = 1.0;
 781                 mLevel[i].gravity = LbmVec(0.0);
 782                 mLevel[i].mprsCells[0] = NULL;
 783                 mLevel[i].mprsCells[1] = NULL;
 784                 mLevel[i].mprsFlags[0] = NULL;
 785                 mLevel[i].mprsFlags[1] = NULL;
 786
 787                 mLevel[i].avgOmega = 0.0;
 788                 mLevel[i].avgOmegaCnt = 0.0;
 789         }
 790
 791         // init sizes
 792         mLevel[mMaxRefine].lSizex = mSizex;
 793         mLevel[mMaxRefine].lSizey = mSizey;
 794         mLevel[mMaxRefine].lSizez = mSizez;
 795         for(int i=mMaxRefine-1; i>=0; i--) {
 796                 mLevel[i].lSizex = mLevel[i+1].lSizex/2;
 797                 mLevel[i].lSizey = mLevel[i+1].lSizey/2;
 798                 mLevel[i].lSizez = mLevel[i+1].lSizez/2;
 799         }
 800
 801         // safety check
 802         if(sizeof(CellFlagType) != CellFlagTypeSize) {
 803                 errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
 804                 return false;
 805         }
 806
 807         double ownMemCheck = 0.0;
 808         mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
 809         mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
 810         mLevel[ mMaxRefine ].lcellfactor = 1.0;
 811         LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
 812
 813 #if COMPRESSGRIDS==0
 814         mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
 815         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
 816         ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
 817 #else // COMPRESSGRIDS==0
 818         LONGINT compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
 819         // D int tmp = ( (rcellSize +compressOffset +4)/(1024*1024) )*4;
 820         // D printf("Debug MEMMMM excee: %d\n", tmp);
 821         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
 822         mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
 823         ownMemCheck += sizeof(LbmFloat) * (rcellSize +compressOffset +4);
 824 #endif // COMPRESSGRIDS==0
 825
 826         if(!mLevel[ mMaxRefine ].mprsCells[1] || !mLevel[ mMaxRefine ].mprsCells[0]) {
 827                 errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (1)! Aborting...",SIMWORLD_INITERROR);
 828                 return false;
 829         }
 830
 831         // +4 for safety ?
 832         mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
 833         mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
 834         ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
 835         if(!mLevel[ mMaxRefine ].mprsFlags[1] || !mLevel[ mMaxRefine ].mprsFlags[0]) {
 836                 errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (2)! Aborting...",SIMWORLD_INITERROR);
 837
 838 #if COMPRESSGRIDS==0
 839                 delete[] mLevel[ mMaxRefine ].mprsCells[0];
 840                 delete[] mLevel[ mMaxRefine ].mprsCells[1];
 841 #else // COMPRESSGRIDS==0
 842                 delete[] mLevel[ mMaxRefine ].mprsCells[1];
 843 #endif // COMPRESSGRIDS==0
 844                 return false;
 845         }
 846
 847         LbmFloat lcfdimFac = 8.0;
 848         if(LBMDIM==2) lcfdimFac = 4.0;
 849         for(int i=mMaxRefine-1; i>=0; i--) {
 850                 mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
 851                 mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
 852                 mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;
 853
 854                 if(LBMDIM==2){ mLevel[i].lSizez = 1; } // 2D
 855                 rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
 856                 mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
 857                 mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
 858                 ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
 859                 mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
 860                 mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
 861                 ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
 862         }
 863
 864         // isosurface memory, use orig res values
 865         if(mFarFieldSize>0.) {
 866                 ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
 867         } else {
 868                 // ignore 3 int slices...
 869                 ownMemCheck += (double)( (              sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
 870         }
 871
 872         // sanity check
 873 #if ELBEEM_PLUGIN!=1
 874         if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
 875                 errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
 876         }
 877 #endif // ELBEEM_PLUGIN!=1
 878
 879         // init sizes for _all_ levels
 880         for(int i=mMaxRefine; i>=0; i--) {
 881                 mLevel[i].lOffsx = mLevel[i].lSizex;
 882                 mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
 883                 mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
 884         mLevel[i].setCurr  = 0;
 885         mLevel[i].setOther = 1;
 886         mLevel[i].lsteps = 0;
 887         mLevel[i].lmass = 0.0;
 888         mLevel[i].lvolume = 0.0;
 889         }
 890
 891         // calc omega, force for all levels
 892         initLevelOmegas();
 893         mMinTimestep = mpParam->getTimestep();
 894         mMaxTimestep = mpParam->getTimestep();
 895
 896         // init isosurf
 897         mpIso->setIsolevel( mIsoValue );
 898 #if LBM_INCLUDE_TESTSOLVERS==1
 899         if(mUseTestdata) {
 900                 mpTest->setMaterialName( mpIso->getMaterialName() );
 901                 delete mpIso;
 902                 mpIso = mpTest;
 903                 if(mpTest->mFarfMode>0) { // 3d off
 904                         mpTest->setIsolevel(-100.0);
 905                 } else {
 906                         mpTest->setIsolevel( mIsoValue );
 907                 }
 908         }
 909 #endif // LBM_INCLUDE_TESTSOLVERS!=1
 910         // approximate feature size with mesh resolution
 911         float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
 912         // smooth vars defined in solver_interface, set by simulation object
 913         // reset for invalid values...
 914         if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
 915         if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
 916         mpIso->setSmoothSurface( mSmoothSurface * featureSize );
 917         mpIso->setSmoothNormals( mSmoothNormals * featureSize );
 918
 919         // init iso weight values mIsoWeightMethod
 920         int wcnt = 0;
 921         float totw = 0.0;
 922         for(int ak=-1;ak<=1;ak++)
 923                 for(int aj=-1;aj<=1;aj++)
 924                         for(int ai=-1;ai<=1;ai++)  {
 925                                 switch(mIsoWeightMethod) {
 926                                 case 1: // light smoothing
 927                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
 928                                         break;
 929                                 case 2: // very light smoothing
 930                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
 931                                         mIsoWeight[wcnt] *= mIsoWeight[wcnt];
 932                                         break;
 933                                 case 3: // no smoothing
 934                                         if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
 935                                         else mIsoWeight[wcnt] = 0.0;
 936                                         break;
 937                                 default: // strong smoothing (=0)
 938                                         mIsoWeight[wcnt] = 1.0;
 939                                         break;
 940                                 }
 941                                 totw += mIsoWeight[wcnt];
 942                                 wcnt++;
 943                         }
 944         for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
 945
 946         LbmVec isostart = vec2L(mvGeoStart);
 947         LbmVec isoend   = vec2L(mvGeoEnd);
 948         int twodOff = 0; // 2d slices
 949         if(LBMDIM==2) {
 950                 LbmFloat sn,se;
 951                 sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
 952                 se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
 953                 isostart[2] = sn;
 954                 isoend[2] = se;
 955                 twodOff = 2;
 956         }
 957         int isosx = mSizex+2;
 958         int isosy = mSizey+2;
 959         int isosz = mSizez+2+twodOff;
 960
 961         // MPT
 962 #if LBM_INCLUDE_TESTSOLVERS==1
 963         //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
 964         if( (mMpNum>0) && (mMpIndex==0) ) {
 965                 //? mpindex==0
 966                 // restore original value for node0
 967                 isosx       = mOrgSizeX + 2;
 968                 isostart[0] = mOrgStartX;
 969                 isoend[0]   = mOrgEndX;
 970         }
 971         errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
 972                         << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
 973                         << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
 974         errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
 975                         << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
 976                         << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
 977 #endif // LBM_INCLUDE_TESTSOLVERS==1
 978
 979         errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
 980         errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
 981                         " grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
 982         mpIso->setStart( vec2G(isostart) );
 983         mpIso->setEnd(   vec2G(isoend) );
 984         LbmVec isodist = isoend-isostart;
 985
 986         int isosubs = mIsoSubdivs;
 987         if(mFarFieldSize>1.) {
 988                 errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
 989                 isosubs = 1;
 990                 mpIso->setUseFulledgeArrays(true);
 991         }
 992         mpIso->setSubdivs(isosubs);
 993
 994         mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
 995
 996         // reset iso field
 997         for(int ak=0;ak<isosz;ak++)
 998                 for(int aj=0;aj<isosy;aj++)
 999                         for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
1000
1001
1002   /* init array (set all invalid first) */
1003         preinitGrids();
1004         for(int lev=0; lev<=mMaxRefine; lev++) {
1005                 FSGR_FORIJK_BOUNDS(lev) {
1006                         RFLAG(lev,i,j,k,0) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
1007                         if(!mAllfluid) {
1008                                 initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0);
1009                         } else {
1010                                 initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0);
1011                         }
1012                 }
1013         }
1014
1015
1016         if(LBMDIM==2) {
1017                 if(mOutputSurfacePreview) {
1018                         errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
1019                         mOutputSurfacePreview = 0; }
1020         }
1021         if((glob_mpactive) && (glob_mpindex>0)) {
1022                 mOutputSurfacePreview = 0;
1023         }
1024
1025 #if LBM_USE_GUI==1
1026         if(mOutputSurfacePreview) {
1027                 errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
1028                 mOutputSurfacePreview = 0; }
1029 #endif // LBM_USE_GUI==1
1030         if(mOutputSurfacePreview) {
1031                 // same as normal one, but use reduced size
1032                 mpPreviewSurface = new IsoSurface( mIsoValue );
1033                 mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
1034                 mpPreviewSurface->setIsolevel( mIsoValue );
1035                 // usually dont display for rendering
1036                 mpPreviewSurface->setVisible( false );
1037
1038                 mpPreviewSurface->setStart( vec2G(isostart) );
1039                 mpPreviewSurface->setEnd(   vec2G(isoend) );
1040                 LbmVec pisodist = isoend-isostart;
1041                 LbmFloat pfac = mPreviewFactor;
1042                 mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
1043                 //mpPreviewSurface->setName( getName() + "preview" );
1044                 mpPreviewSurface->setName( "preview" );
1045
1046                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
1047         }
1048
1049         // init defaults
1050         mAvgNumUsedCells = 0;
1051         mFixMass= 0.0;
1052         return true;
1053 }
1054
1055 /*! init solver arrays */
1056 bool LbmFsgrSolver::initializeSolverGrids() {
1057   /* init boundaries */
1058   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
1059         // init obstacles, and reinit time step size
1060         initGeometryFlags();
1061         mLastSimTime = -1.0;
1062         // TODO check for invalid cells? nitGenericTestCases();
1063
1064         // new - init noslip 1 everywhere...
1065         // half fill boundary cells?
1066
1067         CellFlagType domainBoundType = CFInvalid;
1068         // TODO use normal object types instad...
1069         if(mDomainBound.find(string("free")) != string::npos) {
1070                 domainBoundType = CFBnd | CFBndFreeslip;
1071         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
1072         } else if(mDomainBound.find(string("part")) != string::npos) {
1073                 domainBoundType = CFBnd | CFBndPartslip; // part slip type
1074         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
1075         } else {
1076                 domainBoundType = CFBnd | CFBndNoslip;
1077         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
1078         }
1079
1080         // use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
1081         int domainobj = (int)(mpGiObjects->size());
1082         domainBoundType |= (domainobj<<24);
1083         //for(int i=0; i<(int)(domainobj+0); i++) {
1084                 //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
1085                 //if((*mpGiObjects)[i] == mpIso) { //check...
1086                 //}
1087         //}
1088         //errMsg("GEOIN"," dm "<<(domainBoundType>>24));
1089
1090   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1091     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
1092                         initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL);
1093                         initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL);
1094     }
1095
1096   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1097     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1098                         initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL);
1099                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL);
1100                         // DEBUG BORDER!
1101                         //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
1102     }
1103
1104         if(LBMDIM == 3) {
1105                 // only for 3D
1106                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1107                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
1108                                 initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL);
1109                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL);
1110                         }
1111         }
1112
1113         // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
1114   /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1115     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1116                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
1117     }
1118   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1119     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
1120                         initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL);
1121     }
1122         // */
1123
1124         /*for(int ii=0; ii<(int)po w_change?(2.0,mMaxRefine)-1; ii++) {
1125                 errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
1126                 for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1127                         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1128                                 initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL);  // SYMM!? 2D?
1129                         }
1130                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1131                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
1132                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL);   // SYMM!? 3D?
1133                         }
1134         }
1135         // Symmetry tests */
1136         // vortt
1137 #if LBM_INCLUDE_TESTSOLVERS==1
1138         if(( strstr( this->getName().c_str(), "vorttfluid" ) != NULL) && (LBMDIM==2)) {
1139                 errMsg("VORTT","init");
1140                 int level=mMaxRefine;
1141                 int cx = mLevel[level].lSizex/2;
1142                 int cyo = mLevel[level].lSizey/2;
1143                 int sx = mLevel[level].lSizex/8;
1144                 int sy = mLevel[level].lSizey/8;
1145                 LbmFloat rho = 1.;
1146                 LbmFloat rhomass = 1.;
1147                 LbmFloat uFactor = 0.15;
1148                 LbmFloat vdist = 1.0;
1149
1150                 int cy1=cyo-(int)(vdist*sy);
1151                 int cy2=cyo+(int)(vdist*sy);
1152
1153                 //for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) {
1154                 for(int j=1;j<mLevel[level].lSizey-1;j++)
1155                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1156                                 LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
1157                                 LbmFloat d2 = norm(LbmVec(cx,cy2,0.)-LbmVec(i,j,0));
1158                                 bool in1 = (d1<=(LbmFloat)(sx));
1159                                 bool in2 = (d2<=(LbmFloat)(sx));
1160                                 LbmVec uvec(0.);
1161                           LbmVec v1 = getNormalized( cross( LbmVec(cx,cy1,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
1162                           LbmVec v2 = getNormalized( cross( LbmVec(cx,cy2,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
1163                                 LbmFloat w1=1., w2=1.;
1164                                 if(!in1) w1=(LbmFloat)(sx)/(1.5*d1);
1165                                 if(!in2) w2=(LbmFloat)(sx)/(1.5*d2);
1166                                 if(!in1) w1=0.; if(!in2) w2=0.; // sharp falloff
1167                           uvec += v1*w1;
1168                           uvec += v2*w2;
1169                                 initVelocityCell(level, i,j,0, CFFluid, rho, rhomass, uvec );
1170                                 //errMsg("VORTT","init uvec"<<uvec);
1171                         }
1172
1173         }
1174 #endif // LBM_INCLUDE_TESTSOLVERS==1
1175
1176         //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
1177
1178         // prepare interface cells
1179         initFreeSurfaces();
1180         initStandingFluidGradient();
1181
1182         // perform first step to init initial mass
1183         mInitialMass = 0.0;
1184         int inmCellCnt = 0;
1185         FSGR_FORIJK1(mMaxRefine) {
1186                 if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
1187                         LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0);
1188                         FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
1189                         mInitialMass += fluidRho;
1190                         inmCellCnt ++;
1191                 } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
1192                         mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
1193                         inmCellCnt ++;
1194                 }
1195         }
1196         mCurrentVolume = mCurrentMass = mInitialMass;
1197
1198         ParamVec cspv = mpParam->calculateCellSize();
1199         if(LBMDIM==2) cspv[2] = 1.0;
1200         inmCellCnt = 1;
1201         double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
1202         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
1203         mInitialMass = 0.0; // reset, and use actual value after first step
1204
1205         //mStartSymm = false;
1206 #if ELBEEM_PLUGIN!=1
1207         if((LBMDIM==2)&&(mSizex<200)) {
1208                 if(!checkSymmetry("init")) {
1209                         errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
1210                 } else {
1211                         errMsg("LbmFsgrSolver::initialize","Symmetric init!");
1212                 }
1213         }
1214 #endif // ELBEEM_PLUGIN!=1
1215         return true;
1216 }
1217
1218
1219 /*! prepare actual simulation start, setup viz etc */
1220 bool LbmFsgrSolver::initializeSolverPostinit() {
1221         // coarsen region
1222         myTime_t fsgrtstart = getTime();
1223         for(int lev=mMaxRefine-1; lev>=0; lev--) {
1224                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
1225                 adaptGrid(lev);
1226                 coarseRestrictFromFine(lev);
1227                 adaptGrid(lev);
1228                 coarseRestrictFromFine(lev);
1229         }
1230         markedClearList();
1231         myTime_t fsgrtend = getTime();
1232         if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
1233         mNumFsgrChanges = 0;
1234
1235         for(int l=0; l<cDirNum; l++) {
1236                 LbmFloat area = 0.5 * 0.5 *0.5;
1237                 if(LBMDIM==2) area = 0.5 * 0.5;
1238
1239                 if(dfVecX[l]!=0) area *= 0.5;
1240                 if(dfVecY[l]!=0) area *= 0.5;
1241                 if(dfVecZ[l]!=0) area *= 0.5;
1242                 mFsgrCellArea[l] = area;
1243         } // l
1244
1245         // make sure both sets are ok
1246         // copy from other to curr
1247         for(int lev=0; lev<=mMaxRefine; lev++) {
1248         FSGR_FORIJK_BOUNDS(lev) {
1249                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1250         } } // first copy flags */
1251
1252
1253         // old mpPreviewSurface init
1254         //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
1255         // make sure fill fracs are right for first surface generation
1256         stepMain();
1257
1258         // prepare once...
1259         mpIso->setParticles(mpParticles, mPartDropMassSub);
1260   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
1261         prepareVisualization();
1262         // copy again for stats counting
1263         for(int lev=0; lev<=mMaxRefine; lev++) {
1264         FSGR_FORIJK_BOUNDS(lev) {
1265                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1266         } } // first copy flags */
1267
1268
1269         // now really done...
1270   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
1271   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
1272         mInitDone = 1;
1273
1274 #if LBM_INCLUDE_CONTROL==1
1275         initCpdata();
1276 #endif // LBM_INCLUDE_CONTROL==1
1277
1278 #if LBM_INCLUDE_TESTSOLVERS==1
1279         initTestdata();
1280 #endif // ELBEEM_PLUGIN!=1
1281         // not inited? dont use...
1282         if(mCutoff<0) mCutoff=0;
1283
1284         initParticles();
1285         return true;
1286 }
1287
1288
1289
1290 // macros for mov obj init
1291 #if LBMDIM==2
1292
1293 #define POS2GRID_CHECK(vec,n) \
1294                                 monTotal++;\
1295                                 int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
1296                                 if(k!=0) continue; \
1297                                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
1298                                 if(i<=0) continue; \
1299                                 if(i>=mLevel[level].lSizex-1) continue; \
1300                                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
1301                                 if(j<=0) continue; \
1302                                 if(j>=mLevel[level].lSizey-1) continue;  \
1303
1304 #else // LBMDIM -> 3
1305 #define POS2GRID_CHECK(vec,n) \
1306                                 monTotal++;\
1307                                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
1308                                 if(i<=0) continue; \
1309                                 if(i>=mLevel[level].lSizex-1) continue; \
1310                                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
1311                                 if(j<=0) continue; \
1312                                 if(j>=mLevel[level].lSizey-1) continue; \
1313                                 const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
1314                                 if(k<=0) continue; \
1315                                 if(k>=mLevel[level].lSizez-1) continue; \
1316
1317 #endif // LBMDIM
1318
1319 // calculate object velocity from vert arrays in objvel vec
1320 #define OBJVEL_CALC \
1321                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
1322                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
1323                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
1324                                 if(usqr>maxusqr) { \
1325                                         /* cutoff at maxVelVal */ \
1326                                         for(int jj=0; jj<3; jj++) { \
1327                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
1328                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
1329                                         } \
1330                                 } } \
1331                                 if(ntype&(CFBndFreeslip)) { \
1332                                         const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
1333                                         const LbmVec oldov=objvel; /*DEBUG*/ \
1334                                         objvel = vec2L((*pNormals)[n]) *dp; \
1335                                         /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
1336                                 } \
1337                                 else if(ntype&(CFBndPartslip)) { \
1338                                         const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
1339                                         const LbmVec oldov=objvel; /*DEBUG*/ \
1340                                         /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
1341                                         const LbmFloat partv = mObjectPartslips[OId]; \
1342                                         /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
1343                                         /* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
1344                                         objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
1345                                 }
1346
1347 #define TTT \
1348
1349
1350 /*****************************************************************************/
1351 //! init moving obstacles for next sim step sim
1352 /*****************************************************************************/
1353 void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
1354         myTime_t monstart = getTime();
1355
1356         // movobj init
1357         const int level = mMaxRefine;
1358         const int workSet = mLevel[level].setCurr;
1359         const int otherSet = mLevel[level].setOther;
1360         LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
1361         // for debugging - check targetTime check during DEFAULT STREAM
1362         LbmFloat targetTime = mSimulationTime + mpParam->getTimestep();
1363         if(mLastSimTime == targetTime) {
1364                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
1365                 return;
1366         }
1367         //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
1368         //if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);
1369
1370         const LbmFloat maxVelVal = 0.1666;
1371         const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
1372
1373         LbmFloat rhomass = 0.0;
1374         CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
1375         CellFlagType ntype = CFInvalid;
1376         // WARNING - copied from geo init!
1377         int numobjs = (int)(mpGiObjects->size());
1378         ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1379         // 2d display as rectangles
1380         ntlVec3Gfx iniPos(0.0);
1381         if(LBMDIM==2) {
1382                 dvec[2] = 1.0;
1383                 iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
1384         } else {
1385                 iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
1386         }
1387
1388         if( (int)mObjectMassMovnd.size() < numobjs) {
1389                 for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
1390                         mObjectMassMovnd.push_back(0.);
1391                 }
1392         }
1393
1394         // stats
1395         int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
1396         int nbored;
1397         for(int OId=0; OId<numobjs; OId++) {
1398                 ntlGeometryObject *obj = (*mpGiObjects)[OId];
1399                 bool skip = false;
1400                 if(obj->getGeoInitId() != mLbmInitId) skip=true;
1401                 if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
1402                 if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
1403                 if(skip) continue;
1404                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);
1405
1406                 if( (obj->getGeoInitType()&FGI_ALLBOUNDS) ||
1407                                 (obj->getGeoInitType()&FGI_FLUID) && staticInit ) {
1408
1409                         otype = ntype = CFInvalid;
1410                         switch(obj->getGeoInitType()) {
1411                                 /* case FGI_BNDPART: // old, use noslip for moving part/free objs
1412                                 case FGI_BNDFREE:
1413                                         if(!staticInit) {
1414                                                 errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
1415                                                 otype = ntype = CFBnd|CFBndNoslip;
1416                                         } else {
1417                                                 if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
1418                                                 if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
1419                                         }
1420                                         break;
1421                                         // off */
1422                                 case FGI_BNDPART: rhomass = BND_FILL;
1423                                         otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
1424                                         break;
1425                                 case FGI_BNDFREE: rhomass = BND_FILL;
1426                                         otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
1427                                         break;
1428                                         // off */
1429                                 case FGI_BNDNO:   rhomass = BND_FILL;
1430                                         otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
1431                                         break;
1432                                 case FGI_FLUID:
1433                                         otype = ntype = CFFluid;
1434                                         break;
1435                                 case FGI_MBNDINFLOW:
1436                                         otype = ntype = CFMbndInflow;
1437                                         break;
1438                                 case FGI_MBNDOUTFLOW:
1439                                         otype = ntype = CFMbndOutflow;
1440                                         break;
1441                         }
1442                         int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
1443                         int active =    ((obj->getGeoActive(targetTime)>0.)? 1:0);
1444                         //errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<"  s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
1445                         // skip inactive in/out flows
1446                         if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
1447                         if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;
1448
1449                         // copied from  recalculateObjectSpeeds
1450                         mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
1451                         debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
1452
1453                         //vector<ntlVec3Gfx> tNormals;
1454                         vector<ntlVec3Gfx> *pNormals = NULL;
1455                         mMOINormals.clear();
1456                         if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
1457
1458                         mMOIVertices.clear();
1459                         if(obj->getMeshAnimated()) {
1460                                 // do two full update
1461                                 // TODO tNormals handling!?
1462                                 mMOIVerticesOld.clear();
1463                                 obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
1464                                 monTrafo += mMOIVerticesOld.size();
1465                                 obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
1466                                 monTrafo += mMOIVertices.size();
1467                                 obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
1468                         } else {
1469                                 // only do transform update
1470                                 obj->getMovingPoints(mMOIVertices,pNormals);
1471                                 mMOIVerticesOld = mMOIVertices;
1472                                 // WARNING - assumes mSimulationTime is global!?
1473                                 obj->applyTransformation(targetTime, &mMOIVertices,pNormals, 0, mMOIVertices.size(), false );
1474                                 monTrafo += mMOIVertices.size();
1475
1476                                 // correct flags from last position, but extrapolate
1477                                 // velocity to next timestep
1478                                 obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
1479                                 monTrafo += mMOIVerticesOld.size();
1480                         }
1481
1482                         // object types
1483                         if(ntype&CFBnd){
1484
1485                                 // check if object is moving at all
1486                                 if(obj->getIsAnimated()) {
1487                                         ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
1488                                         // FIXME?
1489                                         if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
1490                                         // get old type - CHECK FIXME , timestep could have changed - cause trouble?
1491                                         ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
1492                                         if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
1493                                 }
1494                                 if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
1495                                 CellFlagType rflagnb[27];
1496                                 LbmFloat massCheck = 0.;
1497                                 int massReinits=0;
1498                                 bool fillCells = (mObjectMassMovnd[OId]<=-1.);
1499                                 LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
1500
1501
1502                                 // first pass - set new obs. cells
1503                                 if(active) {
1504                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1505                                                 //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
1506                                                 POS2GRID_CHECK(mMOIVertices,n);
1507                                                 //{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
1508                                                 if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
1509                                                 monPoints++;
1510
1511                                                 // check mass
1512                                                 if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
1513                                                         FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
1514                                                         massReinits++;
1515                                                 }
1516                                                 else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
1517                                                         massCheck -= QCELL(level, i,j,k, workSet, dMass);
1518                                                         massReinits++;
1519                                                 }
1520
1521                                                 RFLAG(level, i,j,k, workSet) = ntype;
1522                                                 FORDF1 {
1523                                                         //CellFlagType flag = RFLAG_NB(level, i,j,k,workSet,l);
1524                                                         rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
1525                                                         if(rflagnb[l]&(CFFluid|CFInter)) {
1526                                                                 rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
1527                                                                 RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l];
1528                                                         }
1529                                                 }
1530                                                 LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
1531                                                 RAC(dstCell,0) = 0.0;
1532                                                 if(ntype&CFBndMoving) {
1533                                                         OBJVEL_CALC;
1534
1535                                                         // compute fluid acceleration
1536                                                         FORDF1 {
1537                                                                 if(rflagnb[l]&(CFFluid|CFInter)) {
1538                                                                         const LbmFloat ux = dfDvecX[l]*objvel[0];
1539                                                                         const LbmFloat uy = dfDvecY[l]*objvel[1];
1540                                                                         const LbmFloat uz = dfDvecZ[l]*objvel[2];
1541
1542                                                                         LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); //
1543                                                                         if(ntype&(CFBndFreeslip|CFBndPartslip)) {
1544                                                                                 // missing, diag mass checks...
1545                                                                                 //if(l<=LBMDIM*2) factor *= 1.4142;
1546                                                                                 factor *= 2.0; // TODO, optimize
1547                                                                         } else {
1548                                                                                 factor *= 1.2; // TODO, optimize
1549                                                                         }
1550                                                                         factor *= impactCorrFactor; // use manual tweaking channel
1551                                                                         RAC(dstCell,l) = factor;
1552                                                                         massCheck += factor;
1553                                                                 } else {
1554                                                                         RAC(dstCell,l) = 0.;
1555                                                                 }
1556                                                         }
1557
1558 #if NEWDIRVELMOTEST==1
1559                                                         FORDF1 { RAC(dstCell,l) = 0.; }
1560                                                         RAC(dstCell,dMass)  = objvel[0];
1561                                                         RAC(dstCell,dFfrac) = objvel[1];
1562                                                         RAC(dstCell,dC)     = objvel[2];
1563 #endif // NEWDIRVELMOTEST==1
1564                                                 } else {
1565                                                         FORDF1 { RAC(dstCell,l) = 0.0; }
1566                                                 }
1567                                                 RAC(dstCell, dFlux) = targetTime;
1568                                                 //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" dflt="<<RAC(dstCell, dFlux) );
1569                                                 monObsts++;
1570                                         } // points
1571                                 } // bnd, is active?
1572
1573                                 // second pass, remove old ones
1574                                 // warning - initEmptyCell et. al dont overwrite OId or persist flags...
1575                                 if(wasActive) {
1576                                         for(size_t n=0; n<mMOIVerticesOld.size(); n++) {
1577                                                 POS2GRID_CHECK(mMOIVerticesOld,n);
1578                                                 monPoints++;
1579                                                 if((RFLAG(level, i,j,k, workSet) == otype) &&
1580                                                          (QCELL(level, i,j,k, workSet, dFlux) != targetTime)) {
1581                                                         // from mainloop
1582                                                         nbored = 0;
1583                                                         // TODO: optimize for OPT3D==0
1584                                                         FORDF1 {
1585                                                                 //rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
1586                                                                 rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance
1587                                                                 nbored |= rflagnb[l];
1588                                                         }
1589                                                         CellFlagType settype = CFInvalid;
1590                                                         if(nbored&CFFluid) {
1591                                                                 settype = CFInter|CFNoInterpolSrc;
1592                                                                 rhomass = 1.5;
1593                                                                 if(!fillCells) rhomass = 0.;
1594
1595                                                                 OBJVEL_CALC;
1596                                                                 if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
1597
1598                                                                 // new interpolate values
1599                                                                 LbmFloat avgrho = 0.0;
1600                                                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1601                                                                 interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
1602                                                                 initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
1603                                                                 //errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
1604                                                                 massCheck += rhomass;
1605                                                         } else {
1606                                                                 settype = CFEmpty; rhomass = 0.0;
1607                                                                 initEmptyCell(level, i,j,k, settype, 1., rhomass );
1608                                                         }
1609                                                         monFluids++;
1610                                                         massReinits++;
1611                                                 } // flag & simtime
1612                                         }
1613                                 }  // wasactive
1614
1615                                 // only compute mass transfer when init is done
1616                                 if(mInitDone) {
1617                                         errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
1618                                                 " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass );
1619                                         mObjectMassMovnd[OId] += massCheck;
1620                                 }
1621                         } // bnd, active
1622
1623                         else if(ntype&CFFluid){
1624                                 // second static init pass
1625                                 if(staticInit) {
1626                                         //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" verts"<<mMOIVertices.size() ,9);
1627                                         CellFlagType setflflag = CFFluid; //|(OId<<24);
1628                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1629                                                 POS2GRID_CHECK(mMOIVertices,n);
1630                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1631                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1632                                                         //changeFlag(level, i,j,k, workSet, setflflag);
1633                                                         initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
1634                                                 }
1635                                                 //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
1636                                         }
1637                                 } // second static inflow pass
1638                         } // fluid
1639
1640                         else if(ntype&CFMbndInflow){
1641                                 // inflow pass - add new fluid cells
1642                                 // this is slightly different for standing inflows,
1643                                 // as the fluid is forced to the desired velocity inside as
1644                                 // well...
1645                                 const LbmFloat iniRho = 1.0;
1646                                 const LbmVec vel(mObjectSpeeds[OId]);
1647                                 const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5;
1648                                 USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
1649                                 //errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
1650
1651                                 for(size_t n=0; n<mMOIVertices.size(); n++) {
1652                                         POS2GRID_CHECK(mMOIVertices,n);
1653                                         // TODO - also reinit interface cells !?
1654                                         if(RFLAG(level, i,j,k, workSet)!=CFEmpty) {
1655                                                 // test prevent particle gen for inflow inter cells
1656                                                 if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag
1657                                                 continue; }
1658                                         monFluids++;
1659
1660                                         // TODO add OPT3D treatment
1661                                         LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
1662                                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
1663                                         RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
1664                                         RAC(tcel, dFlux) = FLUX_INIT;
1665                                         CellFlagType setFlag = CFInter;
1666                                         changeFlag(level, i,j,k, workSet, setFlag);
1667                                         mInitialMass += iniRho;
1668                                 }
1669                                 // second static init pass
1670                                 if(staticInit) {
1671                                         CellFlagType set2Flag = CFMbndInflow|(OId<<24);
1672                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1673                                                 POS2GRID_CHECK(mMOIVertices,n);
1674                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1675                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1676                                                         forceChangeFlag(level, i,j,k, workSet, set2Flag);
1677                                                 } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
1678                                                         forceChangeFlag(level, i,j,k, workSet,
1679                                                                         (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
1680                                                 }
1681                                         }
1682                                 } // second static inflow pass
1683
1684                         } // inflow
1685
1686                         else if(ntype&CFMbndOutflow){
1687                                 const LbmFloat iniRho = 0.0;
1688                                 for(size_t n=0; n<mMOIVertices.size(); n++) {
1689                                         POS2GRID_CHECK(mMOIVertices,n);
1690                                         // FIXME check fluid/inter cells for non-static!?
1691                                         if(!(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter))) {
1692                                                 if((staticInit)&&(RFLAG(level, i,j,k, workSet)==CFEmpty)) {
1693                                                         forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
1694                                                 continue;
1695                                         }
1696                                         monFluids++;
1697                                         // interface cell - might be double empty? FIXME check?
1698
1699                                         // remove fluid cells, shouldnt be here anyway
1700                                         LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
1701                                         LbmFloat fluidRho = RAC(tcel,0); FORDF1 { fluidRho += RAC(tcel,l); }
1702                                         mInitialMass -= fluidRho;
1703                                         RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
1704                                         RAC(tcel, dFlux) = FLUX_INIT;
1705                                         CellFlagType setFlag = CFInter;
1706                                         changeFlag(level, i,j,k, workSet, setFlag);
1707
1708                                         // same as ifemptied for if below
1709                                         LbmPoint emptyp;
1710                                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1711                                         mListEmpty.push_back( emptyp );
1712                                         //calcCellsEmptied++;
1713                                 } // points
1714                                 // second static init pass
1715                                 if(staticInit) {
1716                                         CellFlagType set2Flag = CFMbndOutflow|(OId<<24);
1717                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1718                                                 POS2GRID_CHECK(mMOIVertices,n);
1719                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1720                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1721                                                         forceChangeFlag(level, i,j,k, workSet, set2Flag);
1722                                                 } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
1723                                                         forceChangeFlag(level, i,j,k, workSet,
1724                                                                         (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
1725                                                 }
1726                                         }
1727                                 } // second static outflow pass
1728                         } // outflow
1729
1730                 } // allbound check
1731         } // OId
1732
1733
1734         /* { // DEBUG check
1735                 int workSet = mLevel[level].setCurr;
1736                 FSGR_FORIJK1(level) {
1737                         if( (RFLAG(level,i,j,k, workSet) & CFBndMoving) ) {
1738                                 if(QCELL(level, i,j,k, workSet, dFlux)!=targetTime) {
1739                                         errMsg("lastt"," old val!? at "<<PRINT_IJK<<" t="<<QCELL(level, i,j,k, workSet, dFlux)<<" target="<<targetTime);
1740                                 }
1741                         }
1742                 }
1743         } // DEBUG */
1744
1745 #undef POS2GRID_CHECK
1746         myTime_t monend = getTime();
1747         if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
1748         mLastSimTime = targetTime;
1749 }
1750
1751
1752 // geoinit position
1753
1754 #define GETPOS(i,j,k)  \
1755         ntlVec3Gfx ggpos = \
1756                 ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
1757                                   iniPos[1]+ dvec[1]*(gfxReal)(j), \
1758                                   iniPos[2]+ dvec[2]*(gfxReal)(k) ); \
1759   if((LBMDIM==2)&&(mInit2dYZ)) { SWAPYZ(ggpos); }
1760
1761 /*****************************************************************************/
1762 /*! perform geometry init (if switched on) */
1763 /*****************************************************************************/
1764 extern int globGeoInitDebug; //solver_interface
1765 bool LbmFsgrSolver::initGeometryFlags() {
1766         int level = mMaxRefine;
1767         myTime_t geotimestart = getTime();
1768         ntlGeometryObject *pObj;
1769         ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1770         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
1771         // WARNING - copied to movobj init!
1772
1773         /* init object velocities, this has always to be called for init */
1774         initGeoTree();
1775         if(mAllfluid) {
1776                 freeGeoTree();
1777                 return true; }
1778
1779         // make sure moving obstacles are inited correctly
1780         // preinit moving obj points
1781         int numobjs = (int)(mpGiObjects->size());
1782         for(int o=0; o<numobjs; o++) {
1783                 ntlGeometryObject *obj = (*mpGiObjects)[o];
1784                 //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
1785                 if(
1786                                 ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
1787                                 (obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
1788                         if(!obj->getMeshAnimated()) {
1789                                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
1790                                 obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
1791                         }
1792                 }
1793         }
1794
1795         // max speed init
1796         ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
1797         ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
1798         mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
1799         LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
1800         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
1801         if(mpParam->getSimulationMaxSpeed() > allowMax) {
1802                 // similar to adaptTimestep();
1803                 LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
1804                 LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
1805                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
1806                 mpParam->setDesiredTimestep( newdt );
1807                 mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
1808                 maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
1809                                       mSimulationTime, mpParam->getTimestep() )) ));
1810                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
1811         }
1812         recalculateObjectSpeeds();
1813         // */
1814
1815         // init obstacles for first time step (requires obj speeds)
1816         initMovingObstacles(true);
1817
1818         /* set interface cells */
1819         ntlVec3Gfx pos,iniPos; // position of current cell
1820         LbmFloat rhomass = 0.0;
1821         CellFlagType ntype = CFInvalid;
1822         int savedNodes = 0;
1823         int OId = -1;
1824         gfxReal distance;
1825
1826         // 2d display as rectangles
1827         if(LBMDIM==2) {
1828                 dvec[2] = 0.0;
1829                 iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
1830                 //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
1831         } else {
1832                 iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
1833         }
1834
1835
1836         // first init boundary conditions
1837         // invalid cells are set to empty afterwards
1838         // TODO use floop macros!?
1839         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1840                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1841                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1842                                 ntype = CFInvalid;
1843
1844                                 GETPOS(i,j,k);
1845                                 const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
1846                                 if(inside) {
1847                                         pObj = (*mpGiObjects)[OId];
1848                                         switch( pObj->getGeoInitType() ){
1849                                         case FGI_MBNDINFLOW:
1850                                           if(! pObj->getIsAnimated() ) {
1851                                                         rhomass = 1.0;
1852                                                         ntype = CFFluid | CFMbndInflow;
1853                                                 } else {
1854                                                         ntype = CFInvalid;
1855                                                 }
1856                                                 break;
1857                                         case FGI_MBNDOUTFLOW:
1858                                           if(! pObj->getIsAnimated() ) {
1859                                                         rhomass = 0.0;
1860                                                         ntype = CFEmpty|CFMbndOutflow;
1861                                                 } else {
1862                                                         ntype = CFInvalid;
1863                                                 }
1864                                                 break;
1865                                         case FGI_BNDNO:
1866                                                 rhomass = BND_FILL;
1867                                                 ntype = CFBnd|CFBndNoslip;
1868                                                 break;
1869                                         case FGI_BNDPART:
1870                                                 rhomass = BND_FILL;
1871                                                 ntype = CFBnd|CFBndPartslip; break;
1872                                         case FGI_BNDFREE:
1873                                                 rhomass = BND_FILL;
1874                                                 ntype = CFBnd|CFBndFreeslip; break;
1875                                         default: // warn here?
1876                                                 rhomass = BND_FILL;
1877                                                 ntype = CFBnd|CFBndNoslip; break;
1878                                         }
1879                                 }
1880                                 if(ntype != CFInvalid) {
1881                                         // initDefaultCell
1882                                         if((ntype & CFMbndInflow) || (ntype & CFMbndOutflow) ) { }
1883                                         ntype |= (OId<<24); // NEWTEST2
1884                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1885                                 }
1886
1887                                 // walk along x until hit for following inits
1888                                 if(distance<=-1.0) { distance = 100.0; } // FIXME dangerous
1889                                 if(distance>=0.0) {
1890                                         gfxReal dcnt=dvec[0];
1891                                         while(( dcnt< distance-dvec[0] )&&(i+1<mLevel[level].lSizex-1)) {
1892                                                 dcnt += dvec[0]; i++;
1893                                                 savedNodes++;
1894                                                 if(ntype != CFInvalid) {
1895                                                         // rho,mass,OId are still inited from above
1896                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1897                                                 }
1898                                         }
1899                                 }
1900                                 // */
1901
1902                         }
1903                 }
1904         } // zmax
1905         // */
1906
1907         // now init fluid layer
1908         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1909                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1910                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1911                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
1912                                 ntype = CFInvalid;
1913                                 int inits = 0;
1914                                 GETPOS(i,j,k);
1915                                 const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
1916                                 if(inside) {
1917                                         ntype = CFFluid;
1918                                 }
1919                                 if(ntype != CFInvalid) {
1920                                         // initDefaultCell
1921                                         rhomass = 1.0;
1922                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1923                                         inits++;
1924                                 }
1925
1926                                 // walk along x until hit for following inits
1927                                 if(distance<=-1.0) { distance = 100.0; }
1928                                 if(distance>=0.0) {
1929                                         gfxReal dcnt=dvec[0];
1930                                         while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
1931                                                 dcnt += dvec[0]; i++;
1932                                                 savedNodes++;
1933                                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
1934                                                 if(ntype != CFInvalid) {
1935                                                         // rhomass are still inited from above
1936                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1937                                                         inits++;
1938                                                 }
1939                                         }
1940                                 } // distance>0
1941
1942                         }
1943                 }
1944         } // zmax
1945
1946         // reset invalid to empty again
1947         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1948                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1949                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1950                                 if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) {
1951                                         RFLAG(level, i,j,k, mLevel[level].setOther) =
1952                                         RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty;
1953                                 }
1954                         }
1955                 }
1956         }
1957
1958         freeGeoTree();
1959         myTime_t geotimeend = getTime();
1960         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 );
1961         //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
1962         return true;
1963 }
1964
1965 #undef GETPOS
1966
1967 /*****************************************************************************/
1968 /* init part for all freesurface testcases */
1969 void LbmFsgrSolver::initFreeSurfaces() {
1970         double interfaceFill = 0.45;   // filling level of interface cells
1971         //interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!
1972
1973         // set interface cells
1974         FSGR_FORIJK1(mMaxRefine) {
1975                 if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
1976                         FORDF1 {
1977                                 int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
1978                                 if( ( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr)& CFEmpty ) ) {
1979                                         LbmFloat arho=0., aux=0., auy=0., auz=0.;
1980                                         interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
1981                                         //errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
1982                                         // unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
1983                                         initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
1984                                 }
1985                         }
1986                 }
1987         }
1988
1989         // remove invalid interface cells
1990         FSGR_FORIJK1(mMaxRefine) {
1991                 // remove invalid interface cells
1992                 if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
1993                         int delit = 0;
1994                         int NBs = 0; // neighbor flags or'ed
1995                         int noEmptyNB = 1;
1996
1997                         FORDF1 {
1998                                 if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
1999                                         noEmptyNB = 0;
2000                                 }
2001                                 NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
2002                         }
2003                         // remove cells with no fluid or interface neighbors
2004                         if((NBs & CFFluid)==0) { delit = 1; }
2005                         if((NBs & CFInter)==0) { delit = 1; }
2006                         // remove cells with no empty neighbors
2007                         if(noEmptyNB) { delit = 2; }
2008
2009                         // now we can remove the cell
2010                         if(delit==1) {
2011                                 initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
2012                         }
2013                         if(delit==2) {
2014                                 //initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
2015                                 LbmFloat arho=0., aux=0., auy=0., auz=0.;
2016                                 interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
2017                                 initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
2018                         }
2019                 } // interface
2020         } // */
2021
2022         // another brute force init, make sure the fill values are right...
2023         // and make sure both sets are equal
2024         for(int lev=0; lev<=mMaxRefine; lev++) {
2025         FSGR_FORIJK_BOUNDS(lev) {
2026                 if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) {
2027                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
2028                         continue;
2029                 }
2030                 if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) {
2031                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
2032                         continue;
2033                 }
2034         } }
2035
2036         // ----------------------------------------------------------------------
2037         // smoother surface...
2038         if(mInitSurfaceSmoothing>0) {
2039                 debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing timestep ",10);
2040 #if COMPRESSGRIDS==1
2041                 //errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
2042 #endif // COMPRESSGRIDS==0
2043         }
2044         for(int s=0; s<mInitSurfaceSmoothing; s++) {
2045                 //SGR_FORIJK1(mMaxRefine) {
2046
2047                 int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
2048                 int lev = mMaxRefine;
2049 #if COMPRESSGRIDS==0
2050                 for(int k=kstart;k<kend;++k) {
2051 #else // COMPRESSGRIDS==0
2052                 int kdir = 1; // COMPRT ON
2053                 if(mLevel[lev].setCurr==1) {
2054                         kdir = -1;
2055                         int temp = kend;
2056                         kend = kstart-1;
2057                         kstart = temp-1;
2058                 } // COMPRT
2059                 for(int k=kstart;k!=kend;k+=kdir) {
2060 #endif // COMPRESSGRIDS==0
2061                 for(int j=1;j<mLevel[lev].lSizey-1;++j) {
2062                 for(int i=1;i<mLevel[lev].lSizex-1;++i) {
2063                         if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
2064                                 LbmFloat mass = 0.0;
2065                                 //LbmFloat nbdiv;
2066                                 //FORDF0 {
2067                                 for(int l=0;(l<cDirNum); l++) {
2068                                         int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
2069                                         if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
2070                                                 mass += 1.0;
2071                                         }
2072                                         if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFInter ){
2073                                                 mass += QCELL(lev, ni,nj,nk, mLevel[lev].setCurr, dMass);
2074                                         }
2075                                         //nbdiv+=1.0;
2076                                 }
2077
2078                                 //errMsg(" I ", PRINT_IJK<<" m"<<mass );
2079                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
2080                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
2081                         }
2082                 }}}
2083
2084                 mLevel[lev].setOther = mLevel[lev].setCurr;
2085                 mLevel[lev].setCurr ^= 1;
2086         }
2087         // copy back...?
2088 }
2089
2090 /*****************************************************************************/
2091 /* init part for all freesurface testcases */
2092 void LbmFsgrSolver::initStandingFluidGradient() {
2093         // ----------------------------------------------------------------------
2094         // standing fluid preinit
2095         const int debugStandingPreinit = 0;
2096         int haveStandingFluid = 0;
2097
2098 #define STANDFLAGCHECK(iindex) \
2099                                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
2100                                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
2101                                         if((iindex)>1) { haveStandingFluid=(iindex); } \
2102                                         j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
2103                                         continue; \
2104                                 }
2105         int gravIndex[3] = {0,0,0};
2106         int gravDir[3] = {1,1,1};
2107         int maxGravComp = 1; // by default y
2108         int gravComp1 = 0; // by default y
2109         int gravComp2 = 2; // by default y
2110         if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
2111         if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
2112
2113         int gravIMin[3] = { 0 , 0 , 0 };
2114         int gravIMax[3] = {
2115                 mLevel[mMaxRefine].lSizex + 0,
2116                 mLevel[mMaxRefine].lSizey + 0,
2117                 mLevel[mMaxRefine].lSizez + 0 };
2118         if(LBMDIM==2) gravIMax[2] = 1;
2119
2120         //int gravDir = 1;
2121         if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
2122                 // swap directions
2123                 int i=maxGravComp;
2124                 int tmp = gravIMin[i];
2125                 gravIMin[i] = gravIMax[i] - 1;
2126                 gravIMax[i] = tmp - 1;
2127                 gravDir[i] = -1;
2128         }
2129 #define PRINTGDIRS \
2130         errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
2131         errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
2132         errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] );
2133         // _PRINTGDIRS;
2134
2135         bool gravAbort = false;
2136 #define GRAVLOOP \
2137         gravAbort=false; \
2138         for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
2139                 for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
2140                         for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0])
2141
2142         GRAVLOOP {
2143                 int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2144                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) ||
2145                     ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) ||
2146                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) {
2147                         int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
2148                         if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
2149                         if(fluidHeight>1) {
2150                                 haveStandingFluid = fluidHeight; //gravIndex[maxGravComp];
2151                                 gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
2152                         }
2153                         gravAbort = true; continue;
2154                 }
2155         } // GRAVLOOP
2156         // _PRINTGDIRS;
2157
2158         LbmFloat fluidHeight;
2159         fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
2160 #if LBM_INCLUDE_TESTSOLVERS==1
2161         if(mUseTestdata) mpTest->mFluidHeight = (int)fluidHeight;
2162 #endif // ELBEEM_PLUGIN!=1
2163         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])<<
2164                         " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
2165
2166         if(mDisableStandingFluidInit) {
2167                 debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
2168                 haveStandingFluid=0;
2169         }
2170
2171         // copy flags and init , as no flags will be changed during grav init
2172         // also important for corasening later on
2173         const int lev = mMaxRefine;
2174         CellFlagType nbflag[LBM_DFNUM], nbored;
2175         for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
2176                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2177                         for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2178                                 if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
2179                                         nbored = 0;
2180                                         FORDF1 {
2181                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2182                                                 nbored |= nbflag[l];
2183                                         }
2184                                         if(nbored&CFBnd) {
2185                                                 RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
2186                                         } else {
2187                                                 RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
2188                                         }
2189                                 }
2190                                 RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
2191         } } }
2192
2193         if(haveStandingFluid) {
2194                 int rhoworkSet = mLevel[lev].setCurr;
2195                 myTime_t timestart = getTime(); // FIXME use user time here?
2196
2197                 GRAVLOOP {
2198                         int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2199                         //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
2200                         if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
2201                                         ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){
2202                                 //gravAbort = true;
2203                                 continue;
2204                         }
2205
2206                         LbmFloat rho = 1.0;
2207                         // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
2208                         rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) *
2209                                 (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega);
2210                         if(debugStandingPreinit)
2211                                 if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) {
2212                                         errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK);
2213                                 }
2214
2215                         if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
2216                                         (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
2217                                 FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
2218                                 QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
2219                         }
2220
2221                 } // GRAVLOOP
2222
2223                 debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<<
2224                         (1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8);
2225
2226                 int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
2227                 preinitSteps = (haveStandingFluid>>2); // not much use...?
2228                 //preinitSteps = 0;
2229                 debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
2230                 for(int s=0; s<preinitSteps; s++) {
2231                         // in solver main cpp
2232                         standingFluidPreinit();
2233                 }
2234
2235                 myTime_t timeend = getTime();
2236                 debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
2237 #undef  NBFLAG
2238         }
2239 }
2240
2241
2242
2243 bool LbmFsgrSolver::checkSymmetry(string idstring)
2244 {
2245         bool erro = false;
2246         bool symm = true;
2247         int msgs = 0;
2248         const int maxMsgs = 10;
2249         const bool markCells = false;
2250
2251         //for(int lev=0; lev<=mMaxRefine; lev++) {
2252         { int lev = mMaxRefine;
2253
2254         // no point if not symm.
2255         if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
2256                 // ok
2257         } else {
2258                 return false;
2259         }
2260
2261         for(int s=0; s<2; s++) {
2262         FSGR_FORIJK1(lev) {
2263                 if(i<(mLevel[lev].lSizex/2)) {
2264                         int inb = (mLevel[lev].lSizey-1-i);
2265
2266                         if(lev==mMaxRefine) inb -= 1;           // FSGR_SYMM_T
2267
2268                         if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
2269                                 if(LBMDIM==2) {
2270                                         if(msgs<maxMsgs) { msgs++;
2271                                                 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) );
2272                                         }
2273                                 }
2274                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2275                                 symm = false;
2276                         }
2277                         if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
2278                                 if(LBMDIM==2) {
2279                                         if(msgs<maxMsgs) { msgs++;
2280                                                 //debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
2281                                                 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) );
2282                                         }
2283                                 }
2284                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2285                                 symm = false;
2286                         }
2287
2288                         LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
2289                         FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
2290                         LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
2291                         FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
2292                         if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
2293                                 if(LBMDIM==2) {
2294                                         if(msgs<maxMsgs) { msgs++;
2295                                                 //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
2296                                                 errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
2297                                         }
2298                                 }
2299                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2300                                 symm = false;
2301                         }
2302                 }
2303         } }
2304         } // lev
2305         LbmFloat maxdiv =0.0;
2306         if(erro) {
2307                 errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
2308                 //if(LBMDIM==2) mPanic = true;
2309                 //return false;
2310         } else {
2311                 errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
2312         }
2313         // all ok...
2314         return symm;
2315 }// */
2316
2317
2318 void
2319 LbmFsgrSolver::interpolateCellValues(
2320                 int level,int ei,int ej,int ek,int workSet,
2321                 LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz)
2322 {
2323         LbmFloat avgrho = 0.0;
2324         LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
2325         LbmFloat cellcnt = 0.0;
2326
2327         for(int nbl=1; nbl< cDfNum ; ++nbl) {
2328                 if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
2329                 if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){
2330                                 //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
2331                                  //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) {
2332                         cellcnt += 1.0;
2333                         for(int rl=0; rl< cDfNum ; ++rl) {
2334                                 LbmFloat nbdf =  QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
2335                                 avgux  += (dfDvecX[rl]*nbdf);
2336                                 avguy  += (dfDvecY[rl]*nbdf);
2337                                 avguz  += (dfDvecZ[rl]*nbdf);
2338                                 avgrho += nbdf;
2339                         }
2340                 }
2341         }
2342
2343         if(cellcnt<=0.0) {
2344                 // no nbs? just use eq.
2345                 avgrho = 1.0;
2346                 avgux = avguy = avguz = 0.0;
2347                 //TTT mNumProblems++;
2348 #if ELBEEM_PLUGIN!=1
2349                 //mPanic=1;
2350                 // this can happen for worst case moving obj scenarios...
2351                 errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
2352 #endif // ELBEEM_PLUGIN
2353         } else {
2354                 // init speed
2355                 avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
2356                 avgrho /= cellcnt;
2357         }
2358
2359         retrho = avgrho;
2360         retux = avgux;
2361         retuy = avguy;
2362         retuz = avguz;
2363 } // interpolateCellValues
2364
2365
2366 /******************************************************************************
2367  * instantiation
2368  *****************************************************************************/
2369
2370 //! lbm factory functions
2371 LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }
2372
2373