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