Blender V2.61 - r43446

solver_init.cpp

Go to the documentation of this file.
00001 
00004 /******************************************************************************
00005  *
00006  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
00007  * Copyright 2003-2006 Nils Thuerey
00008  *
00009  * Standard LBM Factory implementation
00010  *
00011  *****************************************************************************/
00012 
00013 
00014 #include "solver_class.h"
00015 #include "solver_relax.h"
00016 // for geo init FGI_ defines
00017 #include "elbeem.h"
00018 
00019 // helper for 2d init
00020 #define SWAPYZ(vec) { \
00021         const LbmFloat tmp = (vec)[2]; \
00022         (vec)[2] = (vec)[1]; (vec)[1] = tmp; }
00023 
00024 
00025 /*****************************************************************************/
00027 
00028 /*****************************************************************************/
00030 #if LBMDIM==3
00031 
00033     const int LbmFsgrSolver::cDimension = 3;
00034 
00035     // Wi factors for collide step 
00036     const LbmFloat LbmFsgrSolver::cCollenZero    = (1.0/3.0);
00037     const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/18.0);
00038     const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
00039 
00041     const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
00042     const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
00043     const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
00044     const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;
00045 
00047     const int    LbmFsgrSolver::cDfNum      = 19;
00049     const int    LbmFsgrSolver::cDirNum     = 27;
00050 
00051     //const string LbmFsgrSolver::dfString[ cDfNum ] = { 
00052     const char* LbmFsgrSolver::dfString[ cDfNum ] = { 
00053         " C", " N"," S"," E"," W"," T"," B",
00054         "NE","NW","SE","SW",
00055         "NT","NB","ST","SB",
00056         "ET","EB","WT","WB"
00057     };
00058 
00059     const int LbmFsgrSolver::dfNorm[ cDfNum ] = { 
00060         cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB, 
00061         cDirNE, cDirNW, cDirSE, cDirSW, 
00062         cDirNT, cDirNB, cDirST, cDirSB, 
00063         cDirET, cDirEB, cDirWT, cDirWB
00064     };
00065 
00066     const int LbmFsgrSolver::dfInv[ cDfNum ] = { 
00067         cDirC,  cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
00068         cDirSW, cDirSE, cDirNW, cDirNE,
00069         cDirSB, cDirST, cDirNB, cDirNT, 
00070         cDirWB, cDirWT, cDirEB, cDirET
00071     };
00072 
00073     const int LbmFsgrSolver::dfRefX[ cDfNum ] = { 
00074         0,  0, 0, 0, 0, 0, 0,
00075         cDirSE, cDirSW, cDirNE, cDirNW,
00076         0, 0, 0, 0, 
00077         cDirEB, cDirET, cDirWB, cDirWT
00078     };
00079 
00080     const int LbmFsgrSolver::dfRefY[ cDfNum ] = { 
00081         0,  0, 0, 0, 0, 0, 0,
00082         cDirNW, cDirNE, cDirSW, cDirSE,
00083         cDirNB, cDirNT, cDirSB, cDirST,
00084         0, 0, 0, 0
00085     };
00086 
00087     const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { 
00088         0,  0, 0, 0, 0, 0, 0,
00089         0, 0, 0, 0, 
00090         cDirST, cDirSB, cDirNT, cDirNB,
00091         cDirWT, cDirWB, cDirET, cDirEB
00092     };
00093 
00094     // Vector Order 3D:
00095     //  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
00096     //  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
00097     //  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
00098     //  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
00099 
00100     const int LbmFsgrSolver::dfVecX[ cDirNum ] = { 
00101         0, 0,0, 1,-1, 0,0,
00102         1,-1,1,-1,
00103         0,0,0,0,
00104         1,1,-1,-1,
00105          1,-1, 1,-1,
00106          1,-1, 1,-1,
00107     };
00108     const int LbmFsgrSolver::dfVecY[ cDirNum ] = { 
00109         0, 1,-1, 0,0,0,0,
00110         1,1,-1,-1,
00111         1,1,-1,-1,
00112         0,0,0,0,
00113          1, 1,-1,-1,
00114          1, 1,-1,-1
00115     };
00116     const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { 
00117         0, 0,0,0,0,1,-1,
00118         0,0,0,0,
00119         1,-1,1,-1,
00120         1,-1,1,-1,
00121          1, 1, 1, 1,
00122         -1,-1,-1,-1
00123     };
00124 
00125     const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
00126         0, 0,0, 1,-1, 0,0,
00127         1,-1,1,-1,
00128         0,0,0,0,
00129         1,1,-1,-1,
00130          1,-1, 1,-1,
00131          1,-1, 1,-1
00132     };
00133     const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
00134         0, 1,-1, 0,0,0,0,
00135         1,1,-1,-1,
00136         1,1,-1,-1,
00137         0,0,0,0,
00138          1, 1,-1,-1,
00139          1, 1,-1,-1
00140     };
00141     const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
00142         0, 0,0,0,0,1,-1,
00143         0,0,0,0,
00144         1,-1,1,-1,
00145         1,-1,1,-1,
00146          1, 1, 1, 1,
00147         -1,-1,-1,-1
00148     };
00149 
00150     /* principal directions */
00151     const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { 
00152         1,-1, 0,0, 0,0
00153     };
00154     const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { 
00155         0,0, 1,-1, 0,0
00156     };
00157     const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { 
00158         0,0, 0,0, 1,-1
00159     };
00160 
00162     LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
00163     LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
00164 
00165 
00166     const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { 
00167         cCollenZero,
00168         cCollenOne, cCollenOne, cCollenOne, 
00169         cCollenOne, cCollenOne, cCollenOne,
00170         cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
00171         cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, 
00172         cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
00173     };
00174 
00175     /* precalculated equilibrium dfs, inited in lbmsolver constructor */
00176     LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
00177 
00178 #else // end LBMDIM==3 , LBMDIM==2
00179 
00180 /*****************************************************************************/
00183 
00184         const int LbmFsgrSolver::cDimension = 2;
00185 
00187         const LbmFloat LbmFsgrSolver::cCollenZero    = (4.0/9.0);
00188         const LbmFloat LbmFsgrSolver::cCollenOne     = (1.0/9.0);
00189         const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
00190 
00192         const LbmFloat LbmFsgrSolver::cMagicNr2    = 1.0005;
00193         const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
00194         const LbmFloat LbmFsgrSolver::cMagicNr     = 1.010001;
00195         const LbmFloat LbmFsgrSolver::cMagicNrNeg  = -0.010001;
00196 
00198         const int LbmFsgrSolver::cDfNum  = 9;
00199         const int LbmFsgrSolver::cDirNum = 9;
00200 
00201     //const string LbmFsgrSolver::dfString[ cDfNum ] = { 
00202     const char* LbmFsgrSolver::dfString[ cDfNum ] = { 
00203         " C", 
00204         " N",   " S", " E", " W",
00205         "NE", "NW", "SE","SW" 
00206     };
00207 
00208     const int LbmFsgrSolver::dfNorm[ cDfNum ] = { 
00209         cDirC, 
00210         cDirN,  cDirS,  cDirE,  cDirW,
00211         cDirNE, cDirNW, cDirSE, cDirSW 
00212     };
00213 
00214     const int LbmFsgrSolver::dfInv[ cDfNum ] = { 
00215         cDirC,  
00216         cDirS,  cDirN,  cDirW,  cDirE,
00217         cDirSW, cDirSE, cDirNW, cDirNE 
00218     };
00219 
00220     const int LbmFsgrSolver::dfRefX[ cDfNum ] = { 
00221         0,  
00222         0,  0,  0,  0,
00223         cDirSE, cDirSW, cDirNE, cDirNW 
00224     };
00225 
00226     const int LbmFsgrSolver::dfRefY[ cDfNum ] = { 
00227         0,  
00228         0,  0,  0,  0,
00229         cDirNW, cDirNE, cDirSW, cDirSE 
00230     };
00231 
00232     const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { 
00233         0,  0, 0, 0, 0,
00234         0, 0, 0, 0
00235     };
00236 
00237     // Vector Order 2D:
00238     // 0  1 2  3  4  5  6 7  8
00239     // 0, 0,0, 1,-1, 1,-1,1,-1 
00240     // 0, 1,-1, 0,0, 1,1,-1,-1 
00241 
00242     const int LbmFsgrSolver::dfVecX[ cDirNum ] = { 
00243         0, 
00244         0,0, 1,-1,
00245         1,-1,1,-1 
00246     };
00247     const int LbmFsgrSolver::dfVecY[ cDirNum ] = { 
00248         0, 
00249         1,-1, 0,0,
00250         1,1,-1,-1 
00251     };
00252     const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { 
00253         0, 0,0,0,0, 0,0,0,0 
00254     };
00255 
00256     const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
00257         0, 
00258         0,0, 1,-1,
00259         1,-1,1,-1 
00260     };
00261     const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
00262         0, 
00263         1,-1, 0,0,
00264         1,1,-1,-1 
00265     };
00266     const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
00267         0, 0,0,0,0, 0,0,0,0 
00268     };
00269 
00270     const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { 
00271         1,-1, 0,0
00272     };
00273     const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { 
00274         0,0, 1,-1
00275     };
00276     const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { 
00277         0,0, 0,0
00278     };
00279 
00280 
00282     LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
00283     LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
00284 
00285 
00286     const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { 
00287         cCollenZero,
00288         cCollenOne, cCollenOne, cCollenOne, cCollenOne, 
00289         cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
00290     };
00291 
00292     /* precalculated equilibrium dfs, inited in lbmsolver constructor */
00293     LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
00294 
00295 // D2Q9 end
00296 #endif  // LBMDIM==2
00297 
00298 
00299 // required globals
00300 extern bool glob_mpactive;
00301 extern int glob_mpnum, glob_mpindex;
00302 
00303 
00304 /******************************************************************************
00305  * Lbm Constructor
00306  *****************************************************************************/
00307 LbmFsgrSolver::LbmFsgrSolver() :
00308     //D(),
00309     mCurrentMass(0.0), mCurrentVolume(0.0),
00310     mNumProblems(0), 
00311     mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
00312     mpPreviewSurface(NULL), 
00313     mTimeAdap(true), mForceTimeStepReduce(false),
00314     mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
00315     mInitSurfaceSmoothing(0), mFsSurfGenSetting(0),
00316     mTimestepReduceLock(0),
00317     mTimeSwitchCounts(0), mTimeMaxvelStepCnt(0),
00318     mSimulationTime(0.0), mLastSimTime(0.0),
00319     mMinTimestep(0.0), mMaxTimestep(0.0),
00320     mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
00321     mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
00322     mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
00323     mIsoWeightMethod(1),
00324     mMaxRefine(1), 
00325     mDfScaleUp(-1.0), mDfScaleDown(-1.0),
00326     mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03
00327     mDebugOmegaRet(0.0),
00328     mLastOmega(1e10), mLastGravity(1e10),
00329     mNumInvIfTotal(0), mNumFsgrChanges(0),
00330     mDisableStandingFluidInit(0),
00331     mInit2dYZ(false),
00332     mForceTadapRefine(-1), mCutoff(-1)
00333 {
00334     mpControl = new LbmControlData();
00335 
00336 #if LBM_INCLUDE_TESTSOLVERS==1
00337     mpTest = new LbmTestdata();
00338     mMpNum = mMpIndex = 0;
00339     mOrgSizeX = 0;
00340     mOrgStartX = 0.;
00341     mOrgEndX = 0.;
00342 #endif // LBM_INCLUDE_TESTSOLVERS!=1
00343     mpIso = new IsoSurface( mIsoValue );
00344 
00345   // init equilibrium dist. func 
00346   LbmFloat rho=1.0;
00347   FORDF0 {
00348         dfEquil[l] = this->getCollideEq( l,rho,  0.0, 0.0, 0.0);
00349   }
00350     dfEquil[dMass] = 1.;
00351     dfEquil[dFfrac] = 1.;
00352     dfEquil[dFlux] = FLUX_INIT;
00353 
00354     // init LES
00355     int odm = 0;
00356     for(int m=0; m<LBMDIM; m++) { 
00357         for(int l=0; l<cDfNum; l++) { 
00358             this->lesCoeffDiag[m][l] = 
00359             this->lesCoeffOffdiag[m][l] = 0.0;
00360         }
00361     }
00362     for(int m=0; m<LBMDIM; m++) { 
00363         for(int n=0; n<LBMDIM; n++) { 
00364             for(int l=1; l<cDfNum; l++) { 
00365                 LbmFloat em;
00366                 switch(m) {
00367                     case 0: em = dfDvecX[l]; break;
00368                     case 1: em = dfDvecY[l]; break;
00369                     case 2: em = dfDvecZ[l]; break;
00370                     default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
00371                 }
00372                 LbmFloat en;
00373                 switch(n) {
00374                     case 0: en = dfDvecX[l]; break;
00375                     case 1: en = dfDvecY[l]; break;
00376                     case 2: en = dfDvecZ[l]; break;
00377                     default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
00378                 }
00379                 const LbmFloat coeff = em*en;
00380                 if(m==n) {
00381                     this->lesCoeffDiag[m][l] = coeff;
00382                 } else {
00383                     if(m>n) {
00384                         this->lesCoeffOffdiag[odm][l] = coeff;
00385                     }
00386                 }
00387             }
00388 
00389             if(m==n) {
00390             } else {
00391                 if(m>n) odm++;
00392             }
00393         }
00394     }
00395 
00396     mDvecNrm[0] = LbmVec(0.0);
00397   FORDF1 {
00398         mDvecNrm[l] = getNormalized( 
00399             LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] ) 
00400             ) * -1.0; 
00401     }
00402 
00403     // calculate gauss weights for restriction
00404     //LbmFloat mGaussw[27];
00405     LbmFloat totGaussw = 0.0;
00406     const LbmFloat alpha = 1.0;
00407     const LbmFloat gw = sqrt(2.0*LBMDIM);
00408 #if ELBEEM_PLUGIN!=1
00409     errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
00410 #endif
00411     for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
00412     //for(int n=0;(n<cDirNum); n++) { 
00413     for(int n=0;(n<cDfNum); n++) { 
00414         const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
00415         LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
00416         mGaussw[n] = w;
00417         totGaussw += w;
00418     }
00419     for(int n=0;(n<cDirNum); n++) { 
00420         mGaussw[n] = mGaussw[n]/totGaussw;
00421     }
00422 
00423 }
00424 
00425 /*****************************************************************************/
00426 /* Destructor */
00427 /*****************************************************************************/
00428 LbmFsgrSolver::~LbmFsgrSolver()
00429 {
00430   if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; }
00431 #if COMPRESSGRIDS==1
00432     delete [] mLevel[mMaxRefine].mprsCells[1];
00433     mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
00434 #endif // COMPRESSGRIDS==1
00435 
00436     for(int i=0; i<=mMaxRefine; i++) {
00437         for(int s=0; s<2; s++) {
00438             if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
00439             if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
00440         }
00441     }
00442     delete mpIso;
00443     if(mpPreviewSurface) delete mpPreviewSurface;
00444     // cleanup done during scene deletion...
00445     
00446     if(mpControl) delete mpControl;
00447 
00448     // always output performance estimate
00449     debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
00450   if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
00451 }
00452 
00453 
00454 
00455 
00456 /******************************************************************************
00457  * initilize variables fom attribute list 
00458  *****************************************************************************/
00459 void LbmFsgrSolver::parseAttrList()
00460 {
00461     LbmSolverInterface::parseStdAttrList();
00462 
00463     string matIso("default");
00464     matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
00465     mpIso->setMaterialName( matIso );
00466     mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
00467     mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
00468     mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false );
00469     mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );
00470 
00471     mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
00472     mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
00473     mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
00474     mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
00475     mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
00476 
00477     // refinement
00478     mMaxRefine = mRefinementDesired;
00479     mMaxRefine  = mpSifAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
00480     if(mMaxRefine<0) mMaxRefine=0;
00481     if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
00482     mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
00483     mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
00484     mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
00485 
00486     // demo mode settings
00487     mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
00488     // FIXME check needed?
00489     mFVArea   = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
00490 
00491     // debugging - skip some time...
00492     double starttimeskip = 0.;
00493     starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
00494     mSimulationTime += starttimeskip;
00495     if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
00496 
00497     mpControl->parseControldataAttrList(mpSifAttrs);
00498 
00499 #if LBM_INCLUDE_TESTSOLVERS==1
00500     mUseTestdata = 0;
00501     mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
00502     mpTest->parseTestdataAttrList(mpSifAttrs);
00503 #ifdef ELBEEM_PLUGIN
00504     mUseTestdata=1; // DEBUG
00505 #endif // ELBEEM_PLUGIN
00506 
00507     mMpNum = mpSifAttrs->readInt("mpnum",  mMpNum ,"LbmFsgrSolver", "mMpNum", false);
00508     mMpIndex = mpSifAttrs->readInt("mpindex",  mMpIndex ,"LbmFsgrSolver", "mMpIndex", false);
00509     if(glob_mpactive) {
00510         // used instead...
00511         mMpNum = glob_mpnum;
00512         mMpIndex = glob_mpindex;
00513     } else {
00514         glob_mpnum = mMpNum;
00515         glob_mpindex = 0;
00516     }
00517     errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<<glob_mpactive<<", "<<glob_mpindex<<"/"<<glob_mpnum);
00518     if(mMpNum>0) {
00519         mUseTestdata=1; // needed in this case...
00520     }
00521 
00522     errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
00523 #else // LBM_INCLUDE_TESTSOLVERS!=1
00524     // not testsolvers, off by default
00525     mUseTestdata = 0;
00526     if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
00527 #endif // LBM_INCLUDE_TESTSOLVERS!=1
00528 
00529     mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
00530     // deprecated!
00531     float mInitialCsmagoCoarse = 0.0;
00532     mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
00533 #if USE_LES==1
00534 #else // USE_LES==1
00535     debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
00536     mInitialCsmago = 0.0;
00537 #endif // USE_LES==1
00538 }
00539 
00540 
00541 /******************************************************************************
00542  * (part of enabling chapter 6 of "Free Surface Flows with Moving and Deforming Objects for LBM")
00543  *****************************************************************************/
00544 void LbmFsgrSolver::setSurfGenSettings(short value)
00545 {
00546     mFsSurfGenSetting = value;
00547 }
00548 
00549 
00550 /******************************************************************************
00551  * Initialize omegas and forces on all levels (for init/timestep change)
00552  *****************************************************************************/
00553 void LbmFsgrSolver::initLevelOmegas()
00554 {
00555     // no explicit settings
00556     mOmega = mpParam->calculateOmega(mSimulationTime);
00557     mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) );
00558     mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused
00559     if(mInit2dYZ) { SWAPYZ(mGravity); }
00560 
00561     // check if last init was ok
00562     LbmFloat gravDelta = norm(mGravity-mLastGravity);
00563     //errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<mOmega<<" - lom:"<<mLastOmega<<" gv:"<<mGravity<<" - "<<mLastGravity<<" , "<<gravDelta  );
00564     if((mOmega == mLastOmega) && (gravDelta<=0.0)) return;
00565 
00566     if(mInitialCsmago<=0.0) {
00567         if(OPT3D==1) {
00568             errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR); 
00569             return;
00570         }
00571     }
00572 
00573     LbmFloat  fineCsmago  = mInitialCsmago;
00574     LbmFloat coarseCsmago = mInitialCsmago;
00575     LbmFloat maxFineCsmago1    = 0.026; 
00576     LbmFloat maxCoarseCsmago1  = 0.029; // try stabilizing
00577     LbmFloat maxFineCsmago2    = 0.028; 
00578     LbmFloat maxCoarseCsmago2  = 0.032; // try stabilizing some more
00579     if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) {
00580         fineCsmago = maxFineCsmago1;
00581         coarseCsmago = maxCoarseCsmago1;
00582     }
00583     if((mMaxRefine>1)&&(mInitialCsmago<maxFineCsmago2)) {
00584         fineCsmago = maxFineCsmago2;
00585         coarseCsmago = maxCoarseCsmago2;
00586     }
00587         
00588 
00589     // use Tau instead of Omega for calculations
00590     { // init base level
00591         int i = mMaxRefine;
00592         mLevel[i].omega    = mOmega;
00593         mLevel[i].timestep = mpParam->getTimestep();
00594         mLevel[i].lcsmago = fineCsmago; //CSMAGO_INITIAL;
00595         mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
00596         mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
00597     }
00598 
00599     // init all sub levels
00600     for(int i=mMaxRefine-1; i>=0; i--) {
00601         //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
00602         double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
00603         nomega                = 1.0/nomega;
00604         mLevel[i].omega       = (LbmFloat)nomega;
00605         mLevel[i].timestep    = 2.0 * mLevel[i+1].timestep;
00606         mLevel[i].lcsmago = coarseCsmago;
00607         mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
00608         mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
00609     }
00610     
00611     // for lbgk
00612     mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
00613     for(int i=mMaxRefine-1; i>=0; i--) {
00614         // should be the same on all levels...
00615         // for lbgk
00616         mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
00617     }
00618 
00619     mLastOmega = mOmega;
00620     mLastGravity = mGravity;
00621     // debug? invalidate old values...
00622     mGravity = -100.0;
00623     mOmega = -100.0;
00624 
00625     for(int i=0; i<=mMaxRefine; i++) {
00626         if(!mSilent) {
00627             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 
00628                     <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
00629                     <<" cmsagp:"<<mLevel[i].lcsmago<<", "
00630                     << " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
00631         } else {
00632             if(!mInitDone) {
00633                 debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
00634                         <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
00635             }
00636         }
00637     }
00638     if(mMaxRefine>0) {
00639         mDfScaleUp   = (mLevel[0  ].timestep/mLevel[0+1].timestep)* (1.0/mLevel[0  ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
00640         mDfScaleDown = (mLevel[0+1].timestep/mLevel[0  ].timestep)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0  ].omega-1.0); // yu
00641     }
00642 }
00643 
00644 
00645 /******************************************************************************
00646  * Init Solver (values should be read from config file)
00647  *****************************************************************************/
00648 
00650 bool LbmFsgrSolver::initializeSolverMemory()
00651 {
00652   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<mInitDone<<" "<<(void*)this,1);
00653 
00654     // init cppf stage
00655     if(mCppfStage>0) {
00656         mSizex *= mCppfStage;
00657         mSizey *= mCppfStage;
00658         mSizez *= mCppfStage;
00659     }
00660     if(mFsSurfGenSetting==-1) {
00661         // all on
00662         mFsSurfGenSetting = 
00663              fssgNormal   | fssgNoNorth  | fssgNoSouth  | fssgNoEast   |
00664              fssgNoWest   | fssgNoTop    | fssgNoBottom | fssgNoObs   ;
00665     }
00666 
00667     // size inits to force cubic cells and mult4 level dimensions
00668     // and make sure we dont allocate too much...
00669     bool memOk = false;
00670     int orgSx = mSizex;
00671     int orgSy = mSizey;
00672     int orgSz = mSizez;
00673     double sizeReduction = 1.0;
00674     double memEstFromFunc = -1.0;
00675     double memEstFine = -1.0;
00676     string memreqStr("");   
00677     bool firstMInit = true;
00678     int minitTries=0;
00679     while(!memOk) {
00680         minitTries++;
00681         initGridSizes( mSizex, mSizey, mSizez,
00682                 mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL);
00683 
00684         // MPT
00685 #if LBM_INCLUDE_TESTSOLVERS==1
00686         if(firstMInit) {
00687             mrSetup();
00688         }
00689 #endif // LBM_INCLUDE_TESTSOLVERS==1
00690         firstMInit=false;
00691 
00692         calculateMemreqEstimate( mSizex, mSizey, mSizez, 
00693                 mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr );
00694         
00695         double memLimit;
00696         string memLimStr("-");
00697         if(sizeof(void*)==4) {
00698             // 32bit system, limit to 2GB
00699             memLimit = 2.0* 1024.0*1024.0*1024.0;
00700             memLimStr = string("2GB");
00701         } else {
00702             // 64bit, just take 16GB as limit for now...
00703             memLimit = 16.0* 1024.0*1024.0*1024.0;
00704             memLimStr = string("16GB");
00705         }
00706 
00707         // restrict max. chunk of 1 mem block to 1GB for windos
00708         bool memBlockAllocProblem = false;
00709         double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
00710         //std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
00711 #ifdef WIN32
00712         double maxWinMemChunk = 1100.*1024.*1024.;
00713         if(sizeof(void *)==4 && memEstFine>maxWinMemChunk) {
00714             memBlockAllocProblem = true;
00715         }
00716 #endif // WIN32
00717 #ifdef __APPLE__
00718         double maxMacMemChunk = 1200.*1024.*1024.;
00719         if(memEstFine> maxMacMemChunk) {
00720             memBlockAllocProblem = true;
00721         }
00722 #endif // Mac
00723         if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
00724             // max memory chunk for 32bit systems 2gig
00725             memBlockAllocProblem = true;
00726         }
00727 
00728         if(memEstFromFunc>memLimit || memBlockAllocProblem) {
00729             sizeReduction *= 0.9;
00730             mSizex = (int)(orgSx * sizeReduction);
00731             mSizey = (int)(orgSy * sizeReduction);
00732             mSizez = (int)(orgSz * sizeReduction);
00733             debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
00734                     //memEstFromFunc<<"/"<<memLimit<<", "<<
00735                     //memEstFine<<"/"<<maxWinMemChunk<<", "<<
00736                     memreqStr<<"/"<<memLimStr<<", "<<
00737                     "retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
00738                     , 3 );
00739         } else {
00740             memOk = true;
00741         } 
00742     }
00743     
00744     mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
00745   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
00746       ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
00747       ", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
00748       ", est. Mem.Req.: "<<memreqStr    ,2);
00749     mpParam->setSize(mSizex, mSizey, mSizez);
00750     if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }
00751 
00752   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
00753         <<"LBM_EPSILON="<<LBM_EPSILON  <<" "
00754         <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG <<" "
00755         <<"OPT3D="<<OPT3D <<" "
00756         <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
00757         <<"MASS_INVALID="<<MASS_INVALID <<" "
00758         <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
00759         <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
00760         <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
00761         <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" " 
00762         <<"USE_LES="<<USE_LES <<" " 
00763         ,10);
00764 
00765     // perform 2D corrections...
00766     if(LBMDIM == 2) mSizez = 1;
00767 
00768     mpParam->setSimulationMaxSpeed(0.0);
00769     if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
00770     mpParam->setTadapLevels( mMaxRefine+1 );
00771 
00772     if(mForceTadapRefine>mMaxRefine) {
00773         mpParam->setTadapLevels( mForceTadapRefine+1 );
00774         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
00775     }
00776 
00777     if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
00778         errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
00779         return false;
00780     }
00781 
00782 
00783     // init vectors
00784     for(int i=0; i<=mMaxRefine; i++) {
00785         mLevel[i].id = i;
00786         mLevel[i].nodeSize = 0.0; 
00787         mLevel[i].simCellSize = 0.0; 
00788         mLevel[i].omega = 0.0; 
00789         mLevel[i].time = 0.0; 
00790         mLevel[i].timestep = 1.0;
00791         mLevel[i].gravity = LbmVec(0.0); 
00792         mLevel[i].mprsCells[0] = NULL;
00793         mLevel[i].mprsCells[1] = NULL;
00794         mLevel[i].mprsFlags[0] = NULL;
00795         mLevel[i].mprsFlags[1] = NULL;
00796 
00797         mLevel[i].avgOmega = 0.0; 
00798         mLevel[i].avgOmegaCnt = 0.0;
00799     }
00800 
00801     // init sizes
00802     mLevel[mMaxRefine].lSizex = mSizex;
00803     mLevel[mMaxRefine].lSizey = mSizey;
00804     mLevel[mMaxRefine].lSizez = mSizez;
00805     for(int i=mMaxRefine-1; i>=0; i--) {
00806         mLevel[i].lSizex = mLevel[i+1].lSizex/2;
00807         mLevel[i].lSizey = mLevel[i+1].lSizey/2;
00808         mLevel[i].lSizez = mLevel[i+1].lSizez/2;
00809     }
00810 
00811     // safety check
00812     if(sizeof(CellFlagType) != CellFlagTypeSize) {
00813         errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
00814         return false;
00815     }
00816 
00817     double ownMemCheck = 0.0;
00818     mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
00819     mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
00820     mLevel[ mMaxRefine ].lcellfactor = 1.0;
00821     LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
00822 
00823 #if COMPRESSGRIDS==0
00824     mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
00825     mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
00826     ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
00827 #else // COMPRESSGRIDS==0
00828     LONGINT compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
00829     // D int tmp = ( (rcellSize +compressOffset +4)/(1024*1024) )*4;
00830     // D printf("Debug MEMMMM excee: %d\n", tmp);
00831     mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
00832     mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
00833     ownMemCheck += sizeof(LbmFloat) * (rcellSize +compressOffset +4);
00834 #endif // COMPRESSGRIDS==0
00835 
00836     if(!mLevel[ mMaxRefine ].mprsCells[1] || !mLevel[ mMaxRefine ].mprsCells[0]) {
00837         errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (1)! Aborting...",SIMWORLD_INITERROR);
00838         return false;
00839     }
00840 
00841     // +4 for safety ?
00842     mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
00843     mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
00844     ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
00845     if(!mLevel[ mMaxRefine ].mprsFlags[1] || !mLevel[ mMaxRefine ].mprsFlags[0]) {
00846         errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (2)! Aborting...",SIMWORLD_INITERROR);
00847 
00848 #if COMPRESSGRIDS==0
00849         delete[] mLevel[ mMaxRefine ].mprsCells[0];
00850         delete[] mLevel[ mMaxRefine ].mprsCells[1];
00851 #else // COMPRESSGRIDS==0
00852         delete[] mLevel[ mMaxRefine ].mprsCells[1];
00853 #endif // COMPRESSGRIDS==0
00854         return false;
00855     }
00856 
00857     LbmFloat lcfdimFac = 8.0;
00858     if(LBMDIM==2) lcfdimFac = 4.0;
00859     for(int i=mMaxRefine-1; i>=0; i--) {
00860         mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
00861         mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
00862         mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;
00863 
00864         if(LBMDIM==2){ mLevel[i].lSizez = 1; } // 2D
00865         rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
00866         mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
00867         mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
00868         ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
00869         mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
00870         mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
00871         ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
00872     }
00873 
00874     // isosurface memory, use orig res values
00875     if(mFarFieldSize>0.) {
00876         ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
00877     } else {
00878         // ignore 3 int slices...
00879         ownMemCheck += (double)( (              sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
00880     }
00881 
00882     // sanity check
00883 #if ELBEEM_PLUGIN!=1
00884     if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
00885         errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
00886     }
00887 #endif // ELBEEM_PLUGIN!=1
00888     
00889     // init sizes for _all_ levels
00890     for(int i=mMaxRefine; i>=0; i--) {
00891         mLevel[i].lOffsx = mLevel[i].lSizex;
00892         mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
00893         mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
00894     mLevel[i].setCurr  = 0;
00895     mLevel[i].setOther = 1;
00896     mLevel[i].lsteps = 0;
00897     mLevel[i].lmass = 0.0;
00898     mLevel[i].lvolume = 0.0;
00899     }
00900 
00901     // calc omega, force for all levels
00902     initLevelOmegas();
00903     mMinTimestep = mpParam->getTimestep();
00904     mMaxTimestep = mpParam->getTimestep();
00905 
00906     // init isosurf
00907     mpIso->setIsolevel( mIsoValue );
00908 #if LBM_INCLUDE_TESTSOLVERS==1
00909     if(mUseTestdata) {
00910         mpTest->setMaterialName( mpIso->getMaterialName() );
00911         delete mpIso;
00912         mpIso = mpTest;
00913         if(mpTest->mFarfMode>0) { // 3d off
00914             mpTest->setIsolevel(-100.0);
00915         } else {
00916             mpTest->setIsolevel( mIsoValue );
00917         }
00918     }
00919 #endif // LBM_INCLUDE_TESTSOLVERS!=1
00920     // approximate feature size with mesh resolution
00921     float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
00922     // smooth vars defined in solver_interface, set by simulation object
00923     // reset for invalid values...
00924     if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
00925     if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
00926     mpIso->setSmoothSurface( mSmoothSurface * featureSize );
00927     mpIso->setSmoothNormals( mSmoothNormals * featureSize );
00928 
00929     // init iso weight values mIsoWeightMethod
00930     int wcnt = 0;
00931     float totw = 0.0;
00932     for(int ak=-1;ak<=1;ak++) 
00933         for(int aj=-1;aj<=1;aj++) 
00934             for(int ai=-1;ai<=1;ai++)  {
00935                 switch(mIsoWeightMethod) {
00936                 case 1: // light smoothing
00937                     mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
00938                     break;
00939                 case 2: // very light smoothing
00940                     mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
00941                     mIsoWeight[wcnt] *= mIsoWeight[wcnt];
00942                     break;
00943                 case 3: // no smoothing
00944                     if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
00945                     else mIsoWeight[wcnt] = 0.0;
00946                     break;
00947                 default: // strong smoothing (=0)
00948                     mIsoWeight[wcnt] = 1.0;
00949                     break;
00950                 }
00951                 totw += mIsoWeight[wcnt];
00952                 wcnt++;
00953             }
00954     for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
00955 
00956     LbmVec isostart = vec2L(mvGeoStart);
00957     LbmVec isoend   = vec2L(mvGeoEnd);
00958     int twodOff = 0; // 2d slices
00959     if(LBMDIM==2) {
00960         LbmFloat sn,se;
00961         sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
00962         se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
00963         isostart[2] = sn;
00964         isoend[2] = se;
00965         twodOff = 2;
00966     }
00967     int isosx = mSizex+2;
00968     int isosy = mSizey+2;
00969     int isosz = mSizez+2+twodOff;
00970 
00971     // MPT
00972 #if LBM_INCLUDE_TESTSOLVERS==1
00973     //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
00974     if( (mMpNum>0) && (mMpIndex==0) ) {
00975         //? mpindex==0
00976         // restore original value for node0
00977         isosx       = mOrgSizeX + 2;
00978         isostart[0] = mOrgStartX;
00979         isoend[0]   = mOrgEndX;
00980     }
00981     errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
00982             << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
00983             << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
00984     errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
00985             << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
00986             << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
00987 #endif // LBM_INCLUDE_TESTSOLVERS==1
00988 
00989     errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
00990     errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
00991             " grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
00992     mpIso->setStart( vec2G(isostart) );
00993     mpIso->setEnd(   vec2G(isoend) );
00994     LbmVec isodist = isoend-isostart;
00995 
00996     int isosubs = mIsoSubdivs;
00997     if(mFarFieldSize>1.) {
00998         errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
00999         isosubs = 1;
01000         mpIso->setUseFulledgeArrays(true);
01001     }
01002     mpIso->setSubdivs(isosubs);
01003 
01004     mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
01005 
01006     // reset iso field
01007     for(int ak=0;ak<isosz;ak++) 
01008         for(int aj=0;aj<isosy;aj++) 
01009             for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
01010 
01011 
01012   /* init array (set all invalid first) */
01013     preinitGrids();
01014     for(int lev=0; lev<=mMaxRefine; lev++) {
01015         FSGR_FORIJK_BOUNDS(lev) {
01016             RFLAG(lev,i,j,k,0) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
01017             if(!mAllfluid) {
01018                 initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
01019             } else {
01020                 initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); 
01021             }
01022         }
01023     }
01024 
01025 
01026     if(LBMDIM==2) {
01027         if(mOutputSurfacePreview) {
01028             errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
01029             mOutputSurfacePreview = 0; }
01030     }
01031     if((glob_mpactive) && (glob_mpindex>0)) {
01032         mOutputSurfacePreview = 0;
01033     }
01034 
01035 #if LBM_USE_GUI==1
01036     if(mOutputSurfacePreview) {
01037         errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
01038         mOutputSurfacePreview = 0; }
01039 #endif // LBM_USE_GUI==1
01040     if(mOutputSurfacePreview) {
01041         // same as normal one, but use reduced size
01042         mpPreviewSurface = new IsoSurface( mIsoValue );
01043         mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
01044         mpPreviewSurface->setIsolevel( mIsoValue );
01045         // usually dont display for rendering
01046         mpPreviewSurface->setVisible( false );
01047 
01048         mpPreviewSurface->setStart( vec2G(isostart) );
01049         mpPreviewSurface->setEnd(   vec2G(isoend) );
01050         LbmVec pisodist = isoend-isostart;
01051         LbmFloat pfac = mPreviewFactor;
01052         mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
01053         //mpPreviewSurface->setName( getName() + "preview" );
01054         mpPreviewSurface->setName( "preview" );
01055     
01056         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
01057     }
01058 
01059     // init defaults
01060     mAvgNumUsedCells = 0;
01061     mFixMass= 0.0;
01062     return true;
01063 }
01064 
01066 bool LbmFsgrSolver::initializeSolverGrids() {
01067   /* init boundaries */
01068   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
01069     // init obstacles, and reinit time step size 
01070     initGeometryFlags();
01071     mLastSimTime = -1.0;
01072     // TODO check for invalid cells? nitGenericTestCases();
01073     
01074     // new - init noslip 1 everywhere...
01075     // half fill boundary cells?
01076 
01077     CellFlagType domainBoundType = CFInvalid;
01078     // TODO use normal object types instad...
01079     if(mDomainBound.find(string("free")) != string::npos) {
01080         domainBoundType = CFBnd | CFBndFreeslip;
01081     debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
01082     } else if(mDomainBound.find(string("part")) != string::npos) {
01083         domainBoundType = CFBnd | CFBndPartslip; // part slip type
01084     debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
01085     } else { 
01086         domainBoundType = CFBnd | CFBndNoslip;
01087     debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
01088     }
01089 
01090     // use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
01091     int domainobj = (int)(mpGiObjects->size());
01092     domainBoundType |= (domainobj<<24);
01093     //for(int i=0; i<(int)(domainobj+0); i++) {
01094         //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
01095         //if((*mpGiObjects)[i] == mpIso) { //check...
01096         //} 
01097     //}
01098     //errMsg("GEOIN"," dm "<<(domainBoundType>>24));
01099 
01100   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
01101     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
01102             initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL); 
01103             initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL); 
01104     }
01105 
01106   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
01107     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
01108             initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL); 
01109             initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL); 
01110             // DEBUG BORDER!
01111             //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
01112     }
01113 
01114     if(LBMDIM == 3) {
01115         // only for 3D
01116         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
01117             for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
01118                 initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL); 
01119                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL); 
01120             }
01121     }
01122 
01123     // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
01124   /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
01125     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
01126             initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
01127     }
01128   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
01129     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
01130             initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL); 
01131     }
01132     // */
01133 
01134     /*for(int ii=0; ii<(int)po w_change?(2.0,mMaxRefine)-1; ii++) {
01135         errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
01136         for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
01137             for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
01138                 initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL);  // SYMM!? 2D?
01139             }
01140         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
01141             for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
01142                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL);   // SYMM!? 3D?
01143             }
01144     }
01145     // Symmetry tests */
01146     // vortt
01147 #if LBM_INCLUDE_TESTSOLVERS==1
01148     if(( strstr( this->getName().c_str(), "vorttfluid" ) != NULL) && (LBMDIM==2)) {
01149         errMsg("VORTT","init");
01150         int level=mMaxRefine;
01151         int cx = mLevel[level].lSizex/2;
01152         int cyo = mLevel[level].lSizey/2;
01153         int sx = mLevel[level].lSizex/8;
01154         int sy = mLevel[level].lSizey/8;
01155         LbmFloat rho = 1.;
01156         LbmFloat rhomass = 1.;
01157         LbmFloat uFactor = 0.15;
01158         LbmFloat vdist = 1.0;
01159 
01160         int cy1=cyo-(int)(vdist*sy);
01161         int cy2=cyo+(int)(vdist*sy);
01162 
01163         //for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) {      
01164         for(int j=1;j<mLevel[level].lSizey-1;j++)
01165             for(int i=1;i<mLevel[level].lSizex-1;i++) {
01166                 LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
01167                 LbmFloat d2 = norm(LbmVec(cx,cy2,0.)-LbmVec(i,j,0));
01168                 bool in1 = (d1<=(LbmFloat)(sx));
01169                 bool in2 = (d2<=(LbmFloat)(sx));
01170                 LbmVec uvec(0.);
01171               LbmVec v1 = getNormalized( cross( LbmVec(cx,cy1,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
01172               LbmVec v2 = getNormalized( cross( LbmVec(cx,cy2,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
01173                 LbmFloat w1=1., w2=1.;
01174                 if(!in1) w1=(LbmFloat)(sx)/(1.5*d1);
01175                 if(!in2) w2=(LbmFloat)(sx)/(1.5*d2);
01176                 if(!in1) w1=0.; if(!in2) w2=0.; // sharp falloff
01177               uvec += v1*w1;
01178               uvec += v2*w2;
01179                 initVelocityCell(level, i,j,0, CFFluid, rho, rhomass, uvec );
01180                 //errMsg("VORTT","init uvec"<<uvec);
01181             }
01182 
01183     }
01184 #endif // LBM_INCLUDE_TESTSOLVERS==1
01185 
01186     //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
01187 
01188     // prepare interface cells
01189     initFreeSurfaces();
01190     initStandingFluidGradient();
01191 
01192     // perform first step to init initial mass
01193     mInitialMass = 0.0;
01194     int inmCellCnt = 0;
01195     FSGR_FORIJK1(mMaxRefine) {
01196         if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
01197             LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); 
01198             FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
01199             mInitialMass += fluidRho;
01200             inmCellCnt ++;
01201         } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
01202             mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
01203             inmCellCnt ++;
01204         }
01205     }
01206     mCurrentVolume = mCurrentMass = mInitialMass;
01207 
01208     ParamVec cspv = mpParam->calculateCellSize();
01209     if(LBMDIM==2) cspv[2] = 1.0;
01210     inmCellCnt = 1;
01211     double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
01212     debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
01213     mInitialMass = 0.0; // reset, and use actual value after first step
01214 
01215     //mStartSymm = false;
01216 #if ELBEEM_PLUGIN!=1
01217     if((LBMDIM==2)&&(mSizex<200)) {
01218         if(!checkSymmetry("init")) {
01219             errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
01220         } else {
01221             errMsg("LbmFsgrSolver::initialize","Symmetric init!");
01222         }
01223     }
01224 #endif // ELBEEM_PLUGIN!=1
01225     return true;
01226 }
01227 
01228 
01230 bool LbmFsgrSolver::initializeSolverPostinit() {
01231     // coarsen region
01232     myTime_t fsgrtstart = getTime(); 
01233     for(int lev=mMaxRefine-1; lev>=0; lev--) {
01234         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
01235         adaptGrid(lev);
01236         coarseRestrictFromFine(lev);
01237         adaptGrid(lev);
01238         coarseRestrictFromFine(lev);
01239     }
01240     markedClearList();
01241     myTime_t fsgrtend = getTime(); 
01242     if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
01243     mNumFsgrChanges = 0;
01244 
01245     for(int l=0; l<cDirNum; l++) { 
01246         LbmFloat area = 0.5 * 0.5 *0.5;
01247         if(LBMDIM==2) area = 0.5 * 0.5;
01248 
01249         if(dfVecX[l]!=0) area *= 0.5;
01250         if(dfVecY[l]!=0) area *= 0.5;
01251         if(dfVecZ[l]!=0) area *= 0.5;
01252         mFsgrCellArea[l] = area;
01253     } // l
01254 
01255     // make sure both sets are ok
01256     // copy from other to curr
01257     for(int lev=0; lev<=mMaxRefine; lev++) {
01258     FSGR_FORIJK_BOUNDS(lev) {
01259         RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
01260     } } // first copy flags */
01261 
01262 
01263     // old mpPreviewSurface init
01264     //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
01265     // make sure fill fracs are right for first surface generation
01266     stepMain();
01267 
01268     // prepare once...
01269     mpIso->setParticles(mpParticles, mPartDropMassSub);
01270   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
01271     prepareVisualization();
01272     // copy again for stats counting
01273     for(int lev=0; lev<=mMaxRefine; lev++) {
01274     FSGR_FORIJK_BOUNDS(lev) {
01275         RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
01276     } } // first copy flags */
01277 
01278 
01279     // now really done...
01280   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
01281   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
01282     mInitDone = 1;
01283 
01284     // init fluid control
01285     initCpdata();
01286 
01287 #if LBM_INCLUDE_TESTSOLVERS==1
01288     initTestdata();
01289 #endif // ELBEEM_PLUGIN!=1
01290     // not inited? dont use...
01291     if(mCutoff<0) mCutoff=0;
01292 
01293     initParticles();
01294     return true;
01295 }
01296 
01297 
01298 
01299 // macros for mov obj init
01300 #if LBMDIM==2
01301 
01302 #define POS2GRID_CHECK(vec,n) \
01303                 monTotal++;\
01304                 int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
01305                 if(k!=0) continue; \
01306                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
01307                 if(i<=0) continue; \
01308                 if(i>=mLevel[level].lSizex-1) continue; \
01309                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
01310                 if(j<=0) continue; \
01311                 if(j>=mLevel[level].lSizey-1) continue;  \
01312 
01313 #else // LBMDIM -> 3
01314 #define POS2GRID_CHECK(vec,n) \
01315                 monTotal++;\
01316                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
01317                 if(i<=0) continue; \
01318                 if(i>=mLevel[level].lSizex-1) continue; \
01319                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
01320                 if(j<=0) continue; \
01321                 if(j>=mLevel[level].lSizey-1) continue; \
01322                 const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
01323                 if(k<=0) continue; \
01324                 if(k>=mLevel[level].lSizez-1) continue; \
01325 
01326 #endif // LBMDIM 
01327 
01328 // calculate object velocity from vert arrays in objvel vec
01329 #define OBJVEL_CALC \
01330                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
01331                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
01332                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
01333                 if(usqr>maxusqr) { \
01334                     /* cutoff at maxVelVal */ \
01335                     for(int jj=0; jj<3; jj++) { \
01336                         if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
01337                         if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
01338                     } \
01339                 } } \
01340                 if(ntype&(CFBndFreeslip)) { \
01341                     const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
01342                     const LbmVec oldov=objvel; /*DEBUG*/ \
01343                     objvel = vec2L((*pNormals)[n]) *dp; \
01344                     /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
01345                 } \
01346                 else if(ntype&(CFBndPartslip)) { \
01347                     const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
01348                     const LbmVec oldov=objvel; /*DEBUG*/ \
01349                     /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
01350                     const LbmFloat partv = mObjectPartslips[OId]; \
01351                     /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
01352                     /* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
01353                     objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
01354                 }
01355 
01356 #define TTT \
01357 
01358 
01359 /*****************************************************************************/
01361 /*****************************************************************************/
01362 void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
01363     myTime_t monstart = getTime();
01364 
01365     // movobj init
01366     const int level = mMaxRefine;
01367     const int workSet = mLevel[level].setCurr;
01368     const int otherSet = mLevel[level].setOther;
01369     LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
01370     // for debugging - check targetTime check during DEFAULT STREAM
01371     LbmFloat targetTime = mSimulationTime + mpParam->getTimestep();
01372     if(mLastSimTime == targetTime) {
01373         debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
01374         return;
01375     }
01376     //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
01377     //if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);
01378 
01379     const LbmFloat maxVelVal = 0.1666;
01380     const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
01381 
01382     LbmFloat rhomass = 0.0;
01383     CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
01384     CellFlagType ntype = CFInvalid;
01385     // WARNING - copied from geo init!
01386     int numobjs = (int)(mpGiObjects->size());
01387     ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
01388     // 2d display as rectangles
01389     ntlVec3Gfx iniPos(0.0);
01390     if(LBMDIM==2) {
01391         dvec[2] = 1.0; 
01392         iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
01393     } else {
01394         iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
01395     }
01396     
01397     if( (int)mObjectMassMovnd.size() < numobjs) {
01398         for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
01399             mObjectMassMovnd.push_back(0.);
01400         }
01401     }
01402     
01403     // stats
01404     int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
01405     int nbored;
01406     for(int OId=0; OId<numobjs; OId++) {
01407         ntlGeometryObject *obj = (*mpGiObjects)[OId];
01408         bool skip = false;
01409         if(obj->getGeoInitId() != mLbmInitId) skip=true;
01410         if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
01411         if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
01412         if(skip) continue;
01413         debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);
01414 
01415         if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || 
01416                 ((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) {
01417 
01418             otype = ntype = CFInvalid;
01419             switch(obj->getGeoInitType()) {
01420                 /* case FGI_BNDPART: // old, use noslip for moving part/free objs
01421                 case FGI_BNDFREE: 
01422                     if(!staticInit) {
01423                         errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
01424                         otype = ntype = CFBnd|CFBndNoslip;
01425                     } else {
01426                         if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
01427                         if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
01428                     }
01429                     break; 
01430                     // off */
01431                 case FGI_BNDPART: rhomass = BND_FILL;
01432                     otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
01433                     break;
01434                 case FGI_BNDFREE: rhomass = BND_FILL;
01435                     otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
01436                     break;
01437                     // off */
01438                 case FGI_BNDNO:   rhomass = BND_FILL;
01439                     otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
01440                     break;
01441                 case FGI_FLUID: 
01442                     otype = ntype = CFFluid; 
01443                     break;
01444                 case FGI_MBNDINFLOW: 
01445                     otype = ntype = CFMbndInflow; 
01446                     break;
01447                 case FGI_MBNDOUTFLOW: 
01448                     otype = ntype = CFMbndOutflow; 
01449                     break;
01450             }
01451             int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
01452             int active =    ((obj->getGeoActive(targetTime)>0.)? 1:0);
01453             //errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<"  s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
01454             // skip inactive in/out flows
01455             if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
01456             if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;
01457 
01458             // copied from  recalculateObjectSpeeds
01459             mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
01460             debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
01461 
01462             //vector<ntlVec3Gfx> tNormals;
01463             vector<ntlVec3Gfx> *pNormals = NULL;
01464             mMOINormals.clear();
01465             if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
01466 
01467             mMOIVertices.clear();
01468             if(obj->getMeshAnimated()) { 
01469                 // do two full update
01470                 // TODO tNormals handling!?
01471                 mMOIVerticesOld.clear();
01472                 obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
01473                 monTrafo += mMOIVerticesOld.size();
01474                 obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
01475                 monTrafo += mMOIVertices.size();
01476                 obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
01477             } else {
01478                 // only do transform update
01479                 obj->getMovingPoints(mMOIVertices,pNormals); // mMOIVertices = mCachedMovPoints
01480                 mMOIVerticesOld = mMOIVertices;
01481                 // WARNING - assumes mSimulationTime is global!?
01482                 obj->applyTransformation(targetTime, &mMOIVertices,pNormals, 0, mMOIVertices.size(), false );
01483                 monTrafo += mMOIVertices.size();
01484 
01485                 // correct flags from last position, but extrapolate
01486                 // velocity to next timestep
01487                 obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
01488                 monTrafo += mMOIVerticesOld.size();
01489             }
01490 
01491             // object types
01492             if(ntype&CFBnd){
01493 
01494                 // check if object is moving at all
01495                 if(obj->getIsAnimated()) {
01496                     ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
01497                     // FIXME?
01498                     if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
01499                     // get old type - CHECK FIXME , timestep could have changed - cause trouble?
01500                     ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
01501                     if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
01502                 }
01503                 if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
01504                 CellFlagType rflagnb[27];
01505                 LbmFloat massCheck = 0.;
01506                 int massReinits=0;              
01507                 bool fillCells = (mObjectMassMovnd[OId]<=-1.);
01508                 LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
01509                 
01510 
01511                 // first pass - set new obs. cells
01512                 if(active) {
01513                     for(size_t n=0; n<mMOIVertices.size(); n++) {
01514                         //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
01515                         POS2GRID_CHECK(mMOIVertices,n);
01516                         //{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
01517                         if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
01518                         monPoints++;
01519                         
01520                         // check mass
01521                         if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
01522                             FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
01523                             massReinits++;
01524                         }
01525                         else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
01526                             massCheck -= QCELL(level, i,j,k, workSet, dMass);
01527                             massReinits++;
01528                         }
01529 
01530                         RFLAG(level, i,j,k, workSet) = ntype;
01531                         FORDF1 {
01532                             //CellFlagType flag = RFLAG_NB(level, i,j,k,workSet,l);
01533                             rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
01534                             if(rflagnb[l]&(CFFluid|CFInter)) {
01535                                 rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
01536                                 RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l]; 
01537                             }
01538                         }
01539                         LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
01540                         RAC(dstCell,0) = 0.0;
01541                         if(ntype&CFBndMoving) {
01542                             OBJVEL_CALC;
01543                             
01544                             // compute fluid acceleration
01545                             FORDF1 {
01546                                 if(rflagnb[l]&(CFFluid|CFInter)) { 
01547                                     const LbmFloat ux = dfDvecX[l]*objvel[0];
01548                                     const LbmFloat uy = dfDvecY[l]*objvel[1];
01549                                     const LbmFloat uz = dfDvecZ[l]*objvel[2];
01550 
01551                                     LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // 
01552                                     if(ntype&(CFBndFreeslip|CFBndPartslip)) {
01553                                         // missing, diag mass checks...
01554                                         //if(l<=LBMDIM*2) factor *= 1.4142;
01555                                         factor *= 2.0; // TODO, optimize
01556                                     } else {
01557                                         factor *= 1.2; // TODO, optimize
01558                                     }
01559                                     factor *= impactCorrFactor; // use manual tweaking channel
01560                                     RAC(dstCell,l) = factor;
01561                                     massCheck += factor;
01562                                 } else {
01563                                     RAC(dstCell,l) = 0.;
01564                                 }
01565                             }
01566 
01567 #if NEWDIRVELMOTEST==1
01568                             FORDF1 { RAC(dstCell,l) = 0.; }
01569                             RAC(dstCell,dMass)  = objvel[0];
01570                             RAC(dstCell,dFfrac) = objvel[1];
01571                             RAC(dstCell,dC)     = objvel[2];
01572 #endif // NEWDIRVELMOTEST==1
01573                         } else {
01574                             FORDF1 { RAC(dstCell,l) = 0.0; }
01575                         }
01576                         RAC(dstCell, dFlux) = targetTime;
01577                         //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" dflt="<<RAC(dstCell, dFlux) );
01578                         monObsts++;
01579                     } // points
01580                 } // bnd, is active?
01581 
01582                 // second pass, remove old ones
01583                 // warning - initEmptyCell et. al dont overwrite OId or persist flags...
01584                 if(wasActive) {
01585                     for(size_t n=0; n<mMOIVerticesOld.size(); n++) {
01586                         POS2GRID_CHECK(mMOIVerticesOld,n);
01587                         monPoints++;
01588                         if((RFLAG(level, i,j,k, workSet) == otype) &&
01589                              (QCELL(level, i,j,k, workSet, dFlux) != targetTime)) {
01590                             // from mainloop
01591                             nbored = 0;
01592                             // TODO: optimize for OPT3D==0
01593                             FORDF1 {
01594                                 //rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
01595                                 rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance
01596                                 nbored |= rflagnb[l];
01597                             } 
01598                             CellFlagType settype = CFInvalid;
01599                             if(nbored&CFFluid) {
01600                                 settype = CFInter|CFNoInterpolSrc; 
01601                                 rhomass = 1.5;
01602                                 if(!fillCells) rhomass = 0.;
01603 
01604                                 OBJVEL_CALC;
01605                                 if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
01606 
01607                                 // new interpolate values
01608                                 LbmFloat avgrho = 0.0;
01609                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
01610                                 interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
01611                                 initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
01612                                 //errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
01613                                 massCheck += rhomass;
01614                             } else {
01615                                 settype = CFEmpty; rhomass = 0.0;
01616                                 initEmptyCell(level, i,j,k, settype, 1., rhomass );
01617                             }
01618                             monFluids++;
01619                             massReinits++;
01620                         } // flag & simtime
01621                     }
01622                 }  // wasactive
01623 
01624                 // only compute mass transfer when init is done
01625                 if(mInitDone) {
01626                     errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
01627                         " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass ); 
01628                     mObjectMassMovnd[OId] += massCheck;
01629                 }
01630             } // bnd, active
01631 
01632             else if(ntype&CFFluid){
01633                 // second static init pass
01634                 if(staticInit) {
01635                     //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" verts"<<mMOIVertices.size() ,9);
01636                     CellFlagType setflflag = CFFluid; //|(OId<<24);
01637                     for(size_t n=0; n<mMOIVertices.size(); n++) {
01638                         POS2GRID_CHECK(mMOIVertices,n);
01639                         if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
01640                         if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
01641                             //changeFlag(level, i,j,k, workSet, setflflag);
01642                             initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
01643                         } 
01644                         //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
01645                     }
01646                 } // second static inflow pass
01647             } // fluid
01648 
01649             else if(ntype&CFMbndInflow){
01650                 // inflow pass - add new fluid cells
01651                 // this is slightly different for standing inflows,
01652                 // as the fluid is forced to the desired velocity inside as 
01653                 // well...
01654                 const LbmFloat iniRho = 1.0;
01655                 const LbmVec vel(mObjectSpeeds[OId]);
01656                 const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5;
01657                 USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
01658                 //errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
01659                 
01660                 for(size_t n=0; n<mMOIVertices.size(); n++) {
01661                     POS2GRID_CHECK(mMOIVertices,n);
01662                     // TODO - also reinit interface cells !?
01663                     if(RFLAG(level, i,j,k, workSet)!=CFEmpty) { 
01664                         // test prevent particle gen for inflow inter cells
01665                         if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag
01666                         continue; }
01667                     monFluids++;
01668 
01669                     // TODO add OPT3D treatment
01670                     LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
01671                     FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
01672                     RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
01673                     RAC(tcel, dFlux) = FLUX_INIT;
01674                     CellFlagType setFlag = CFInter;
01675                     changeFlag(level, i,j,k, workSet, setFlag);
01676                     mInitialMass += iniRho;
01677                 }
01678                 // second static init pass
01679                 if(staticInit) {
01680                     CellFlagType set2Flag = CFMbndInflow|(OId<<24);
01681                     for(size_t n=0; n<mMOIVertices.size(); n++) {
01682                         POS2GRID_CHECK(mMOIVertices,n);
01683                         if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
01684                         if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
01685                             forceChangeFlag(level, i,j,k, workSet, set2Flag);
01686                         } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
01687                             forceChangeFlag(level, i,j,k, workSet, 
01688                                     (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
01689                         }
01690                     }
01691                 } // second static inflow pass
01692 
01693             } // inflow
01694 
01695             else if(ntype&CFMbndOutflow){
01696                 const LbmFloat iniRho = 0.0;
01697                 for(size_t n=0; n<mMOIVertices.size(); n++) {
01698                     POS2GRID_CHECK(mMOIVertices,n);
01699                     // FIXME check fluid/inter cells for non-static!?
01700                     if(!(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter))) {
01701                         if((staticInit)&&(RFLAG(level, i,j,k, workSet)==CFEmpty)) {
01702                             forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
01703                         continue;
01704                     }
01705                     monFluids++;
01706                     // interface cell - might be double empty? FIXME check?
01707 
01708                     // remove fluid cells, shouldnt be here anyway
01709                     LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
01710                     LbmFloat fluidRho = RAC(tcel,0); FORDF1 { fluidRho += RAC(tcel,l); }
01711                     mInitialMass -= fluidRho;
01712                     RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
01713                     RAC(tcel, dFlux) = FLUX_INIT;
01714                     CellFlagType setFlag = CFInter;
01715                     changeFlag(level, i,j,k, workSet, setFlag);
01716 
01717                     // same as ifemptied for if below
01718                     LbmPoint emptyp;
01719                     emptyp.x = i; emptyp.y = j; emptyp.z = k;
01720                     mListEmpty.push_back( emptyp );
01721                     //calcCellsEmptied++;
01722                 } // points
01723                 // second static init pass
01724                 if(staticInit) {
01725                     CellFlagType set2Flag = CFMbndOutflow|(OId<<24);
01726                     for(size_t n=0; n<mMOIVertices.size(); n++) {
01727                         POS2GRID_CHECK(mMOIVertices,n);
01728                         if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
01729                         if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
01730                             forceChangeFlag(level, i,j,k, workSet, set2Flag);
01731                         } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
01732                             forceChangeFlag(level, i,j,k, workSet, 
01733                                     (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
01734                         }
01735                     }
01736                 } // second static outflow pass
01737             } // outflow
01738 
01739         } // allbound check
01740     } // OId
01741 
01742 
01743     /* { // DEBUG check
01744         int workSet = mLevel[level].setCurr;
01745         FSGR_FORIJK1(level) {
01746             if( (RFLAG(level,i,j,k, workSet) & CFBndMoving) ) {
01747                 if(QCELL(level, i,j,k, workSet, dFlux)!=targetTime) {
01748                     errMsg("lastt"," old val!? at "<<PRINT_IJK<<" t="<<QCELL(level, i,j,k, workSet, dFlux)<<" target="<<targetTime);
01749                 }
01750             }
01751         }
01752     } // DEBUG */
01753     
01754 #undef POS2GRID_CHECK
01755     myTime_t monend = getTime();
01756     if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
01757     mLastSimTime = targetTime;
01758 }
01759 
01760 
01761 // geoinit position
01762 
01763 #define GETPOS(i,j,k)  \
01764     ntlVec3Gfx ggpos = \
01765         ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
01766                       iniPos[1]+ dvec[1]*(gfxReal)(j), \
01767                       iniPos[2]+ dvec[2]*(gfxReal)(k) ); \
01768   if((LBMDIM==2)&&(mInit2dYZ)) { SWAPYZ(ggpos); }
01769 
01770 /*****************************************************************************/
01772 /*****************************************************************************/
01773 extern int globGeoInitDebug; //solver_interface
01774 bool LbmFsgrSolver::initGeometryFlags() {
01775     int level = mMaxRefine;
01776     myTime_t geotimestart = getTime(); 
01777     ntlGeometryObject *pObj;
01778     ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
01779     debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
01780     // WARNING - copied to movobj init!
01781 
01782     /* init object velocities, this has always to be called for init */
01783     initGeoTree();
01784     if(mAllfluid) { 
01785         freeGeoTree();
01786         return true; }
01787 
01788     // make sure moving obstacles are inited correctly
01789     // preinit moving obj points
01790     int numobjs = (int)(mpGiObjects->size());
01791     for(int o=0; o<numobjs; o++) {
01792         ntlGeometryObject *obj = (*mpGiObjects)[o];
01793         //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
01794         if(
01795                 ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
01796                 (obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
01797             if(!obj->getMeshAnimated()) { 
01798                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
01799                 obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
01800             }
01801         }
01802     }
01803 
01804     // max speed init
01805     ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
01806     ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
01807     mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
01808     LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
01809     debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
01810     if(mpParam->getSimulationMaxSpeed() > allowMax) {
01811         // similar to adaptTimestep();
01812         LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
01813         LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
01814         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
01815         mpParam->setDesiredTimestep( newdt );
01816         mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
01817         maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
01818                               mSimulationTime, mpParam->getTimestep() )) ));
01819         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
01820     }
01821     recalculateObjectSpeeds();
01822     // */
01823 
01824     // init obstacles for first time step (requires obj speeds)
01825     initMovingObstacles(true);
01826 
01827     /* set interface cells */
01828     ntlVec3Gfx pos,iniPos; // position of current cell
01829     LbmFloat rhomass = 0.0;
01830     CellFlagType ntype = CFInvalid;
01831     int savedNodes = 0;
01832     int OId = -1;
01833     gfxReal distance;
01834 
01835     // 2d display as rectangles
01836     if(LBMDIM==2) {
01837         dvec[2] = 0.0; 
01838         iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
01839         //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
01840     } else {
01841         iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
01842     }
01843 
01844 
01845     // first init boundary conditions
01846     // invalid cells are set to empty afterwards
01847     // TODO use floop macros!?
01848     for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
01849         for(int j=1;j<mLevel[level].lSizey-1;j++) {
01850             for(int i=1;i<mLevel[level].lSizex-1;i++) {
01851                 ntype = CFInvalid;
01852                 
01853                 GETPOS(i,j,k);
01854                 const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
01855                 if(inside) {
01856                     pObj = (*mpGiObjects)[OId];
01857                     switch( pObj->getGeoInitType() ){
01858                     case FGI_MBNDINFLOW:  
01859                       if(! pObj->getIsAnimated() ) {
01860                             rhomass = 1.0;
01861                             ntype = CFFluid | CFMbndInflow;
01862                         } else {
01863                             ntype = CFInvalid;
01864                         }
01865                         break;
01866                     case FGI_MBNDOUTFLOW: 
01867                       if(! pObj->getIsAnimated() ) {
01868                             rhomass = 0.0;
01869                             ntype = CFEmpty|CFMbndOutflow; 
01870                         } else {
01871                             ntype = CFInvalid;
01872                         }
01873                         break;
01874                     case FGI_BNDNO: 
01875                         rhomass = BND_FILL;
01876                         ntype = CFBnd|CFBndNoslip; 
01877                         break;
01878                     case FGI_BNDPART: 
01879                         rhomass = BND_FILL;
01880                         ntype = CFBnd|CFBndPartslip; break;
01881                     case FGI_BNDFREE: 
01882                         rhomass = BND_FILL;
01883                         ntype = CFBnd|CFBndFreeslip; break;
01884                     default: // warn here?
01885                         rhomass = BND_FILL;
01886                         ntype = CFBnd|CFBndNoslip; break;
01887                     }
01888                 }
01889                 if(ntype != CFInvalid) {
01890                     // initDefaultCell
01891                     if((ntype & CFMbndInflow) || (ntype & CFMbndOutflow) ) { }
01892                     ntype |= (OId<<24); // NEWTEST2
01893                     initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
01894                 }
01895 
01896                 // walk along x until hit for following inits
01897                 if(distance<=-1.0) { distance = 100.0; } // FIXME dangerous
01898                 if(distance>=0.0) {
01899                     gfxReal dcnt=dvec[0];
01900                     while(( dcnt< distance-dvec[0] )&&(i+1<mLevel[level].lSizex-1)) {
01901                         dcnt += dvec[0]; i++;
01902                         savedNodes++;
01903                         if(ntype != CFInvalid) {
01904                             // rho,mass,OId are still inited from above
01905                             initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
01906                         }
01907                     }
01908                 } 
01909                 // */
01910 
01911             } 
01912         } 
01913     } // zmax
01914     // */
01915 
01916     // now init fluid layer
01917     for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
01918         for(int j=1;j<mLevel[level].lSizey-1;j++) {
01919             for(int i=1;i<mLevel[level].lSizex-1;i++) {
01920                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
01921                 ntype = CFInvalid;
01922                 int inits = 0;
01923                 GETPOS(i,j,k);
01924                 const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
01925                 if(inside) {
01926                     ntype = CFFluid;
01927                 }
01928                 if(ntype != CFInvalid) {
01929                     // initDefaultCell
01930                     rhomass = 1.0;
01931                     initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
01932                     inits++;
01933                 }
01934 
01935                 // walk along x until hit for following inits
01936                 if(distance<=-1.0) { distance = 100.0; }
01937                 if(distance>=0.0) {
01938                     gfxReal dcnt=dvec[0];
01939                     while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
01940                         dcnt += dvec[0]; i++;
01941                         savedNodes++;
01942                         if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
01943                         if(ntype != CFInvalid) {
01944                             // rhomass are still inited from above
01945                             initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
01946                             inits++;
01947                         }
01948                     }
01949                 } // distance>0
01950                 
01951             } 
01952         } 
01953     } // zmax
01954 
01955     // reset invalid to empty again
01956     for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
01957         for(int j=1;j<mLevel[level].lSizey-1;j++) {
01958             for(int i=1;i<mLevel[level].lSizex-1;i++) {
01959                 if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) {
01960                     RFLAG(level, i,j,k, mLevel[level].setOther) = 
01961                     RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty;
01962                 }
01963             }
01964         }
01965     }
01966 
01967     freeGeoTree();
01968     myTime_t geotimeend = getTime(); 
01969     debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 ); 
01970     //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
01971     return true;
01972 }
01973 
01974 #undef GETPOS
01975 
01976 /*****************************************************************************/
01977 /* init part for all freesurface testcases */
01978 void LbmFsgrSolver::initFreeSurfaces() {
01979     double interfaceFill = 0.45;   // filling level of interface cells
01980     //interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!
01981 
01982     // set interface cells 
01983     FSGR_FORIJK1(mMaxRefine) {
01984         if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
01985             FORDF1 {
01986                 int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
01987                 if( ( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr)& CFEmpty ) ) {
01988                     LbmFloat arho=0., aux=0., auy=0., auz=0.;
01989                     interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
01990                     //errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
01991                     // unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
01992                     initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
01993                 }
01994             }
01995         }
01996     }
01997 
01998     // remove invalid interface cells 
01999     FSGR_FORIJK1(mMaxRefine) {
02000         // remove invalid interface cells 
02001         if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
02002             int delit = 0;
02003             int NBs = 0; // neighbor flags or'ed 
02004             int noEmptyNB = 1;
02005 
02006             FORDF1 {
02007                 if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
02008                     noEmptyNB = 0;
02009                 }
02010                 NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
02011             }
02012             // remove cells with no fluid or interface neighbors
02013             if((NBs & CFFluid)==0) { delit = 1; }
02014             if((NBs & CFInter)==0) { delit = 1; }
02015             // remove cells with no empty neighbors
02016             if(noEmptyNB) { delit = 2; }
02017 
02018             // now we can remove the cell 
02019             if(delit==1) {
02020                 initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
02021             }
02022             if(delit==2) {
02023                 //initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
02024                 LbmFloat arho=0., aux=0., auy=0., auz=0.;
02025                 interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
02026                 initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
02027             }
02028         } // interface 
02029     } // */
02030 
02031     // another brute force init, make sure the fill values are right...
02032     // and make sure both sets are equal
02033     for(int lev=0; lev<=mMaxRefine; lev++) {
02034     FSGR_FORIJK_BOUNDS(lev) {
02035         if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { 
02036             QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
02037             continue;
02038         }
02039         if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { 
02040             QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
02041             continue;
02042         }
02043     } }
02044 
02045     // ----------------------------------------------------------------------
02046     // smoother surface...
02047     if(mInitSurfaceSmoothing>0) {
02048         debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing timestep ",10);
02049 #if COMPRESSGRIDS==1
02050         //errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
02051 #endif // COMPRESSGRIDS==0
02052     }
02053     for(int s=0; s<mInitSurfaceSmoothing; s++) {
02054         //SGR_FORIJK1(mMaxRefine) {
02055 
02056         int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
02057         int lev = mMaxRefine;
02058 #if COMPRESSGRIDS==0
02059         for(int k=kstart;k<kend;++k) {
02060 #else // COMPRESSGRIDS==0
02061         int kdir = 1; // COMPRT ON
02062         if(mLevel[lev].setCurr==1) {
02063             kdir = -1;
02064             int temp = kend;
02065             kend = kstart-1;
02066             kstart = temp-1;
02067         } // COMPRT
02068         for(int k=kstart;k!=kend;k+=kdir) {
02069 #endif // COMPRESSGRIDS==0
02070         for(int j=1;j<mLevel[lev].lSizey-1;++j) {
02071         for(int i=1;i<mLevel[lev].lSizex-1;++i) {
02072             if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
02073                 LbmFloat mass = 0.0;
02074                 //LbmFloat nbdiv;
02075                 //FORDF0 {
02076                 for(int l=0;(l<cDirNum); l++) { 
02077                     int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
02078                     if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
02079                         mass += 1.0;
02080                     }
02081                     if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFInter ){
02082                         mass += QCELL(lev, ni,nj,nk, mLevel[lev].setCurr, dMass);
02083                     }
02084                     //nbdiv+=1.0;
02085                 }
02086 
02087                 //errMsg(" I ", PRINT_IJK<<" m"<<mass );
02088                 QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
02089                 QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
02090             }
02091         }}}
02092 
02093         mLevel[lev].setOther = mLevel[lev].setCurr;
02094         mLevel[lev].setCurr ^= 1;
02095     }
02096     // copy back...?
02097 }
02098 
02099 /*****************************************************************************/
02100 /* init part for all freesurface testcases */
02101 void LbmFsgrSolver::initStandingFluidGradient() {
02102     // ----------------------------------------------------------------------
02103     // standing fluid preinit
02104     const int debugStandingPreinit = 0;
02105     int haveStandingFluid = 0;
02106 
02107 #define STANDFLAGCHECK(iindex) \
02108                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
02109                         ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
02110                     if((iindex)>1) { haveStandingFluid=(iindex); } \
02111                     j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
02112                     continue; \
02113                 } 
02114     int gravIndex[3] = {0,0,0};
02115     int gravDir[3] = {1,1,1};
02116     int maxGravComp = 1; // by default y
02117     int gravComp1 = 0; // by default y
02118     int gravComp2 = 2; // by default y
02119     if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
02120     if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
02121 
02122     int gravIMin[3] = { 0 , 0 , 0 };
02123     int gravIMax[3] = {
02124         mLevel[mMaxRefine].lSizex + 0,
02125         mLevel[mMaxRefine].lSizey + 0,
02126         mLevel[mMaxRefine].lSizez + 0 };
02127     if(LBMDIM==2) gravIMax[2] = 1;
02128 
02129     //int gravDir = 1;
02130     if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
02131         // swap directions
02132         int i=maxGravComp;
02133         int tmp = gravIMin[i];
02134         gravIMin[i] = gravIMax[i] - 1;
02135         gravIMax[i] = tmp - 1;
02136         gravDir[i] = -1;
02137     }
02138 #define PRINTGDIRS \
02139     errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
02140     errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
02141     errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); 
02142     // _PRINTGDIRS;
02143 
02144     bool gravAbort = false;
02145 #define GRAVLOOP \
02146     gravAbort=false; \
02147     for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
02148         for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
02149             for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0]) 
02150 
02151     GRAVLOOP {
02152         int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
02153         if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || 
02154             ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) || 
02155                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) {  
02156             int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
02157             if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
02158             if(fluidHeight>1) {
02159                 haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; 
02160                 gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
02161             }
02162             gravAbort = true; continue; 
02163         } 
02164     } // GRAVLOOP
02165     // _PRINTGDIRS;
02166 
02167     LbmFloat fluidHeight;
02168     fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
02169 #if LBM_INCLUDE_TESTSOLVERS==1
02170     if(mUseTestdata) mpTest->mFluidHeight = (int)fluidHeight;
02171 #endif // ELBEEM_PLUGIN!=1
02172     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])<<
02173             " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
02174                 
02175     if(mDisableStandingFluidInit) {
02176         debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
02177         haveStandingFluid=0;
02178     }
02179 
02180     // copy flags and init , as no flags will be changed during grav init
02181     // also important for corasening later on
02182     const int lev = mMaxRefine;
02183     CellFlagType nbflag[LBM_DFNUM], nbored; 
02184     for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
02185         for(int j=0;j<mLevel[lev].lSizey-0;++j) {
02186             for(int i=0;i<mLevel[lev].lSizex-0;++i) {
02187                 if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
02188                     nbored = 0;
02189                     FORDF1 {
02190                         nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
02191                         nbored |= nbflag[l];
02192                     } 
02193                     if(nbored&CFBnd) {
02194                         RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
02195                     } else {
02196                         RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
02197                     }
02198                 }
02199                 RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
02200     } } }
02201 
02202     if(haveStandingFluid) {
02203         int rhoworkSet = mLevel[lev].setCurr;
02204         myTime_t timestart = getTime(); // FIXME use user time here?
02205 
02206         GRAVLOOP {
02207             int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
02208             //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
02209             if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
02210                     ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ 
02211                 //gravAbort = true; 
02212                 continue;
02213             }
02214 
02215             LbmFloat rho = 1.0;
02216             // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
02217             rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * 
02218                 (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
02219             if(debugStandingPreinit) 
02220                 if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { 
02221                     errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); 
02222                 }
02223 
02224             if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
02225                     (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
02226                 FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
02227                 QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
02228             }
02229 
02230         } // GRAVLOOP
02231 
02232         debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<<
02233             (1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8);
02234         
02235         int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
02236         preinitSteps = (haveStandingFluid>>2); // not much use...?
02237         //preinitSteps = 0;
02238         debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
02239         for(int s=0; s<preinitSteps; s++) {
02240             // in solver main cpp
02241             standingFluidPreinit();
02242         }
02243 
02244         myTime_t timeend = getTime();
02245         debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
02246 #undef  NBFLAG
02247     }
02248 }
02249 
02250 
02251 
02252 bool LbmFsgrSolver::checkSymmetry(string idstring)
02253 {
02254     bool erro = false;
02255     bool symm = true;
02256     int msgs = 0;
02257     const int maxMsgs = 10;
02258     const bool markCells = false;
02259 
02260     //for(int lev=0; lev<=mMaxRefine; lev++) {
02261     { int lev = mMaxRefine;
02262 
02263     // no point if not symm.
02264     if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
02265         // ok
02266     } else {
02267         return false;
02268     }
02269 
02270     for(int s=0; s<2; s++) {
02271     FSGR_FORIJK1(lev) {
02272         if(i<(mLevel[lev].lSizex/2)) {
02273             int inb = (mLevel[lev].lSizey-1-i); 
02274 
02275             if(lev==mMaxRefine) inb -= 1;       // FSGR_SYMM_T
02276 
02277             if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
02278                 if(LBMDIM==2) {
02279                     if(msgs<maxMsgs) { msgs++;
02280                         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) );
02281                     }
02282                 }
02283                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
02284                 symm = false;
02285             }
02286             if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
02287                 if(LBMDIM==2) {
02288                     if(msgs<maxMsgs) { msgs++;
02289                         //debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
02290                         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) );
02291                     }
02292                 }
02293                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
02294                 symm = false;
02295             }
02296 
02297             LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
02298             FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
02299             LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
02300             FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
02301             if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
02302                 if(LBMDIM==2) {
02303                     if(msgs<maxMsgs) { msgs++;
02304                         //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
02305                         errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
02306                     }
02307                 }
02308                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
02309                 symm = false;
02310             }
02311         }
02312     } }
02313     } // lev
02314     LbmFloat maxdiv =0.0;
02315     if(erro) {
02316         errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
02317         //if(LBMDIM==2) mPanic = true; 
02318         //return false;
02319     } else {
02320         errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
02321     }
02322     // all ok...
02323     return symm;
02324 }// */
02325 
02326 
02327 void 
02328 LbmFsgrSolver::interpolateCellValues(
02329         int level,int ei,int ej,int ek,int workSet, 
02330         LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz) 
02331 {
02332     LbmFloat avgrho = 0.0;
02333     LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
02334     LbmFloat cellcnt = 0.0;
02335 
02336     for(int nbl=1; nbl< cDfNum ; ++nbl) {
02337         if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
02338         if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){
02339                 //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
02340                  //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) { 
02341             cellcnt += 1.0;
02342             for(int rl=0; rl< cDfNum ; ++rl) { 
02343                 LbmFloat nbdf =  QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
02344                 avgux  += (dfDvecX[rl]*nbdf); 
02345                 avguy  += (dfDvecY[rl]*nbdf);  
02346                 avguz  += (dfDvecZ[rl]*nbdf);  
02347                 avgrho += nbdf;
02348             }
02349         }
02350     }
02351 
02352     if(cellcnt<=0.0) {
02353         // no nbs? just use eq.
02354         avgrho = 1.0;
02355         avgux = avguy = avguz = 0.0;
02356         //TTT mNumProblems++;
02357 #if ELBEEM_PLUGIN!=1
02358         //mPanic=1; 
02359         // this can happen for worst case moving obj scenarios...
02360         errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
02361 #endif // ELBEEM_PLUGIN
02362     } else {
02363         // init speed
02364         avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
02365         avgrho /= cellcnt;
02366     }
02367 
02368     retrho = avgrho;
02369     retux = avgux;
02370     retuy = avguy;
02371     retuz = avguz;
02372 } // interpolateCellValues
02373 
02374 
02375 /******************************************************************************
02376  * instantiation
02377  *****************************************************************************/
02378 
02380 LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }
02381 
02382