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