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