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