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