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