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