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