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