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