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