Bugfix [#34749] Fluid domain > 10GB crashes Blender - out of memory
[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                 double memLimit;
695                 string memLimStr("-");
696                 if(sizeof(void*)==4) {
697                         // 32bit system, limit to 2GB
698                         memLimit = 2.0* 1024.0*1024.0*1024.0;
699                         memLimStr = string("2GB");
700                 } else {
701                         // 64bit, just take 16GB as limit for now...
702                         memLimit = 16.0* 1024.0*1024.0*1024.0;
703                         memLimStr = string("16GB");
704                 }
705
706                 // restrict max. chunk of 1 mem block to 1GB for windos
707                 bool memBlockAllocProblem = false;
708                 double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
709                 //std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
710 #ifdef WIN32
711                 double maxWinMemChunk = 1100.*1024.*1024.;
712                 if(sizeof(void *)==4 && memEstFine>maxWinMemChunk) {
713                         memBlockAllocProblem = true;
714                 }
715 #endif // WIN32
716 #ifdef __APPLE__
717                 double maxMacMemChunk = 1200.*1024.*1024.;
718                 if(memEstFine> maxMacMemChunk) {
719                         memBlockAllocProblem = true;
720                 }
721 #endif // Mac
722                 if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
723                         // max memory chunk for 32bit systems 2gig
724                         memBlockAllocProblem = true;
725                 }
726
727                 if(memEstFromFunc>memLimit || memBlockAllocProblem) {
728                         sizeReduction *= 0.9;
729                         mSizex = (int)(orgSx * sizeReduction);
730                         mSizey = (int)(orgSy * sizeReduction);
731                         mSizez = (int)(orgSz * sizeReduction);
732                         debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
733                                         //memEstFromFunc<<"/"<<memLimit<<", "<<
734                                         //memEstFine<<"/"<<maxWinMemChunk<<", "<<
735                                         memreqStr<<"/"<<memLimStr<<", "<<
736                                         "retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
737                                         , 3 );
738                 } else {
739                         memOk = true;
740                 } 
741         }
742         
743         mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
744   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
745           ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
746           ", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
747           ", est. Mem.Req.: "<<memreqStr        ,2);
748         mpParam->setSize(mSizex, mSizey, mSizez);
749         if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }
750
751   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
752                 <<"LBM_EPSILON="<<LBM_EPSILON  <<" "
753                 <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG <<" "
754                 <<"OPT3D="<<OPT3D <<" "
755                 <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
756                 <<"MASS_INVALID="<<MASS_INVALID <<" "
757                 <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
758                 <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
759                 <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
760                 <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" " 
761                 <<"USE_LES="<<USE_LES <<" " 
762                 ,10);
763
764         // perform 2D corrections...
765         if(LBMDIM == 2) mSizez = 1;
766
767         mpParam->setSimulationMaxSpeed(0.0);
768         if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
769         mpParam->setTadapLevels( mMaxRefine+1 );
770
771         if(mForceTadapRefine>mMaxRefine) {
772                 mpParam->setTadapLevels( mForceTadapRefine+1 );
773                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
774         }
775
776         if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
777                 errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
778                 return false;
779         }
780
781
782         // init vectors
783         for(int i=0; i<=mMaxRefine; i++) {
784                 mLevel[i].id = i;
785                 mLevel[i].nodeSize = 0.0; 
786                 mLevel[i].simCellSize = 0.0; 
787                 mLevel[i].omega = 0.0; 
788                 mLevel[i].time = 0.0; 
789                 mLevel[i].timestep = 1.0;
790                 mLevel[i].gravity = LbmVec(0.0); 
791                 mLevel[i].mprsCells[0] = NULL;
792                 mLevel[i].mprsCells[1] = NULL;
793                 mLevel[i].mprsFlags[0] = NULL;
794                 mLevel[i].mprsFlags[1] = NULL;
795
796                 mLevel[i].avgOmega = 0.0; 
797                 mLevel[i].avgOmegaCnt = 0.0;
798         }
799
800         // init sizes
801         mLevel[mMaxRefine].lSizex = mSizex;
802         mLevel[mMaxRefine].lSizey = mSizey;
803         mLevel[mMaxRefine].lSizez = mSizez;
804         for(int i=mMaxRefine-1; i>=0; i--) {
805                 mLevel[i].lSizex = mLevel[i+1].lSizex/2;
806                 mLevel[i].lSizey = mLevel[i+1].lSizey/2;
807                 mLevel[i].lSizez = mLevel[i+1].lSizez/2;
808         }
809
810         // safety check
811         if(sizeof(CellFlagType) != CellFlagTypeSize) {
812                 errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
813                 return false;
814         }
815
816         double ownMemCheck = 0.0;
817         mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
818         mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
819         mLevel[ mMaxRefine ].lcellfactor = 1.0;
820         LONGINT rcellSize = (LONGINT)((LONGINT)((LONGINT)mLevel[mMaxRefine].lSizex*(LONGINT)mLevel[mMaxRefine].lSizey*(LONGINT)mLevel[mMaxRefine].lSizez) * (LONGINT)dTotalNum);
821
822 #if COMPRESSGRIDS==0
823         mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
824         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
825         ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
826 #else // COMPRESSGRIDS==0
827         LONGINT compressOffset = (LONGINT)((LONGINT)mLevel[mMaxRefine].lSizex * (LONGINT)mLevel[mMaxRefine].lSizey * (LONGINT)dTotalNum * 2);
828         // LONGINT tmp = ( (rcellSize +compressOffset +4)/(1024*1024) )*sizeof(LbmFloat);
829         // printf("Debug MEMMMM excee: %I64d, %I64d, %I64d, %d, %d\n", tmp, compressOffset, rcellSize, mLevel[mMaxRefine].lSizex, mLevel[mMaxRefine].lSizey );
830         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
831         mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
832         ownMemCheck += sizeof(LbmFloat) * (rcellSize +compressOffset +4);
833 #endif // COMPRESSGRIDS==0
834
835         if(!mLevel[ mMaxRefine ].mprsCells[1] || !mLevel[ mMaxRefine ].mprsCells[0]) {
836                 errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (1)! Aborting...",SIMWORLD_INITERROR);
837                 return false;
838         }
839
840         // +4 for safety ?
841         mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
842         mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
843         ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
844         if(!mLevel[ mMaxRefine ].mprsFlags[1] || !mLevel[ mMaxRefine ].mprsFlags[0]) {
845                 errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (2)! Aborting...",SIMWORLD_INITERROR);
846
847 #if COMPRESSGRIDS==0
848                 delete[] mLevel[ mMaxRefine ].mprsCells[0];
849                 delete[] mLevel[ mMaxRefine ].mprsCells[1];
850 #else // COMPRESSGRIDS==0
851                 delete[] mLevel[ mMaxRefine ].mprsCells[1];
852 #endif // COMPRESSGRIDS==0
853                 return false;
854         }
855
856         LbmFloat lcfdimFac = 8.0;
857         if(LBMDIM==2) lcfdimFac = 4.0;
858         for(int i=mMaxRefine-1; i>=0; i--) {
859                 mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
860                 mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
861                 mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;
862
863                 if(LBMDIM==2){ mLevel[i].lSizez = 1; } // 2D
864                 rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
865                 mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
866                 mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
867                 ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
868                 mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
869                 mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
870                 ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
871         }
872
873         // isosurface memory, use orig res values
874         if(mFarFieldSize>0.) {
875                 ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
876         } else {
877                 // ignore 3 int slices...
878                 ownMemCheck += (double)( (              sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
879         }
880
881         // sanity check
882 #if ELBEEM_PLUGIN!=1
883         if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
884                 errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
885         }
886 #endif // ELBEEM_PLUGIN!=1
887         
888         // init sizes for _all_ levels
889         for(int i=mMaxRefine; i>=0; i--) {
890                 mLevel[i].lOffsx = mLevel[i].lSizex;
891                 mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
892                 mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
893         mLevel[i].setCurr  = 0;
894         mLevel[i].setOther = 1;
895         mLevel[i].lsteps = 0;
896         mLevel[i].lmass = 0.0;
897         mLevel[i].lvolume = 0.0;
898         }
899
900         // calc omega, force for all levels
901         initLevelOmegas();
902         mMinTimestep = mpParam->getTimestep();
903         mMaxTimestep = mpParam->getTimestep();
904
905         // init isosurf
906         mpIso->setIsolevel( mIsoValue );
907 #if LBM_INCLUDE_TESTSOLVERS==1
908         if(mUseTestdata) {
909                 mpTest->setMaterialName( mpIso->getMaterialName() );
910                 delete mpIso;
911                 mpIso = mpTest;
912                 if(mpTest->mFarfMode>0) { // 3d off
913                         mpTest->setIsolevel(-100.0);
914                 } else {
915                         mpTest->setIsolevel( mIsoValue );
916                 }
917         }
918 #endif // LBM_INCLUDE_TESTSOLVERS!=1
919         // approximate feature size with mesh resolution
920         float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
921         // smooth vars defined in solver_interface, set by simulation object
922         // reset for invalid values...
923         if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
924         if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
925         mpIso->setSmoothSurface( mSmoothSurface * featureSize );
926         mpIso->setSmoothNormals( mSmoothNormals * featureSize );
927
928         // init iso weight values mIsoWeightMethod
929         int wcnt = 0;
930         float totw = 0.0;
931         for(int ak=-1;ak<=1;ak++) 
932                 for(int aj=-1;aj<=1;aj++) 
933                         for(int ai=-1;ai<=1;ai++)  {
934                                 switch(mIsoWeightMethod) {
935                                 case 1: // light smoothing
936                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
937                                         break;
938                                 case 2: // very light smoothing
939                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
940                                         mIsoWeight[wcnt] *= mIsoWeight[wcnt];
941                                         break;
942                                 case 3: // no smoothing
943                                         if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
944                                         else mIsoWeight[wcnt] = 0.0;
945                                         break;
946                                 default: // strong smoothing (=0)
947                                         mIsoWeight[wcnt] = 1.0;
948                                         break;
949                                 }
950                                 totw += mIsoWeight[wcnt];
951                                 wcnt++;
952                         }
953         for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
954
955         LbmVec isostart = vec2L(mvGeoStart);
956         LbmVec isoend   = vec2L(mvGeoEnd);
957         int twodOff = 0; // 2d slices
958         if(LBMDIM==2) {
959                 LbmFloat sn,se;
960                 sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
961                 se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
962                 isostart[2] = sn;
963                 isoend[2] = se;
964                 twodOff = 2;
965         }
966         int isosx = mSizex+2;
967         int isosy = mSizey+2;
968         int isosz = mSizez+2+twodOff;
969
970         // MPT
971 #if LBM_INCLUDE_TESTSOLVERS==1
972         //if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
973         if( (mMpNum>0) && (mMpIndex==0) ) {
974                 //? mpindex==0
975                 // restore original value for node0
976                 isosx       = mOrgSizeX + 2;
977                 isostart[0] = mOrgStartX;
978                 isoend[0]   = mOrgEndX;
979         }
980         errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
981                         << PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
982                         << PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
983         errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
984                         << PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
985                         << PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
986 #endif // LBM_INCLUDE_TESTSOLVERS==1
987
988         errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
989         errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
990                         " grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
991         mpIso->setStart( vec2G(isostart) );
992         mpIso->setEnd(   vec2G(isoend) );
993         LbmVec isodist = isoend-isostart;
994
995         int isosubs = mIsoSubdivs;
996         if(mFarFieldSize>1.) {
997                 errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
998                 isosubs = 1;
999                 mpIso->setUseFulledgeArrays(true);
1000         }
1001         mpIso->setSubdivs(isosubs);
1002
1003         mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
1004
1005         // reset iso field
1006         for(int ak=0;ak<isosz;ak++) 
1007                 for(int aj=0;aj<isosy;aj++) 
1008                         for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
1009
1010
1011   /* init array (set all invalid first) */
1012         preinitGrids();
1013         for(int lev=0; lev<=mMaxRefine; lev++) {
1014                 FSGR_FORIJK_BOUNDS(lev) {
1015                         RFLAG(lev,i,j,k,0) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
1016                         if(!mAllfluid) {
1017                                 initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
1018                         } else {
1019                                 initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); 
1020                         }
1021                 }
1022         }
1023
1024
1025         if(LBMDIM==2) {
1026                 if(mOutputSurfacePreview) {
1027                         errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
1028                         mOutputSurfacePreview = 0; }
1029         }
1030         if((glob_mpactive) && (glob_mpindex>0)) {
1031                 mOutputSurfacePreview = 0;
1032         }
1033
1034 #if LBM_USE_GUI==1
1035         if(mOutputSurfacePreview) {
1036                 errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
1037                 mOutputSurfacePreview = 0; }
1038 #endif // LBM_USE_GUI==1
1039         if(mOutputSurfacePreview) {
1040                 // same as normal one, but use reduced size
1041                 mpPreviewSurface = new IsoSurface( mIsoValue );
1042                 mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
1043                 mpPreviewSurface->setIsolevel( mIsoValue );
1044                 // usually dont display for rendering
1045                 mpPreviewSurface->setVisible( false );
1046
1047                 mpPreviewSurface->setStart( vec2G(isostart) );
1048                 mpPreviewSurface->setEnd(   vec2G(isoend) );
1049                 LbmVec pisodist = isoend-isostart;
1050                 LbmFloat pfac = mPreviewFactor;
1051                 mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
1052                 //mpPreviewSurface->setName( getName() + "preview" );
1053                 mpPreviewSurface->setName( "preview" );
1054         
1055                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
1056         }
1057
1058         // init defaults
1059         mAvgNumUsedCells = 0;
1060         mFixMass= 0.0;
1061         return true;
1062 }
1063
1064 /*! init solver arrays */
1065 bool LbmFsgrSolver::initializeSolverGrids() {
1066   /* init boundaries */
1067   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
1068         // init obstacles, and reinit time step size 
1069         initGeometryFlags();
1070         mLastSimTime = -1.0;
1071         // TODO check for invalid cells? nitGenericTestCases();
1072         
1073         // new - init noslip 1 everywhere...
1074         // half fill boundary cells?
1075
1076         CellFlagType domainBoundType = CFInvalid;
1077         // TODO use normal object types instad...
1078         if(mDomainBound.find(string("free")) != string::npos) {
1079                 domainBoundType = CFBnd | CFBndFreeslip;
1080         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
1081         } else if(mDomainBound.find(string("part")) != string::npos) {
1082                 domainBoundType = CFBnd | CFBndPartslip; // part slip type
1083         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
1084         } else { 
1085                 domainBoundType = CFBnd | CFBndNoslip;
1086         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
1087         }
1088
1089         // use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
1090         int domainobj = (int)(mpGiObjects->size());
1091         domainBoundType |= (domainobj<<24);
1092         //for(int i=0; i<(int)(domainobj+0); i++) {
1093                 //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
1094                 //if((*mpGiObjects)[i] == mpIso) { //check...
1095                 //} 
1096         //}
1097         //errMsg("GEOIN"," dm "<<(domainBoundType>>24));
1098
1099   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1100     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1101                         initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL); 
1102                         initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL); 
1103     }
1104
1105   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1106     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1107                         initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL); 
1108                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL); 
1109                         // DEBUG BORDER!
1110                         //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
1111     }
1112
1113         if(LBMDIM == 3) {
1114                 // only for 3D
1115                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1116                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1117                                 initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL); 
1118                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL); 
1119                         }
1120         }
1121
1122         // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
1123   /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1124     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1125                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); 
1126     }
1127   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1128     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1129                         initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL); 
1130     }
1131         // */
1132
1133         /*for(int ii=0; ii<(int)po w_change?(2.0,mMaxRefine)-1; ii++) {
1134                 errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
1135                 for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1136                         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1137                                 initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL);  // SYMM!? 2D?
1138                         }
1139                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1140                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1141                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL);   // SYMM!? 3D?
1142                         }
1143         }
1144         // Symmetry tests */
1145         // vortt
1146 #if LBM_INCLUDE_TESTSOLVERS==1
1147         if(( strstr( this->getName().c_str(), "vorttfluid" ) != NULL) && (LBMDIM==2)) {
1148                 errMsg("VORTT","init");
1149                 int level=mMaxRefine;
1150                 int cx = mLevel[level].lSizex/2;
1151                 int cyo = mLevel[level].lSizey/2;
1152                 int sx = mLevel[level].lSizex/8;
1153                 int sy = mLevel[level].lSizey/8;
1154                 LbmFloat rho = 1.;
1155                 LbmFloat rhomass = 1.;
1156                 LbmFloat uFactor = 0.15;
1157                 LbmFloat vdist = 1.0;
1158
1159                 int cy1=cyo-(int)(vdist*sy);
1160                 int cy2=cyo+(int)(vdist*sy);
1161
1162                 //for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) {      
1163                 for(int j=1;j<mLevel[level].lSizey-1;j++)
1164                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1165                                 LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
1166                                 LbmFloat d2 = norm(LbmVec(cx,cy2,0.)-LbmVec(i,j,0));
1167                                 bool in1 = (d1<=(LbmFloat)(sx));
1168                                 bool in2 = (d2<=(LbmFloat)(sx));
1169                                 LbmVec uvec(0.);
1170                           LbmVec v1 = getNormalized( cross( LbmVec(cx,cy1,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
1171                           LbmVec v2 = getNormalized( cross( LbmVec(cx,cy2,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )*  uFactor;
1172                                 LbmFloat w1=1., w2=1.;
1173                                 if(!in1) w1=(LbmFloat)(sx)/(1.5*d1);
1174                                 if(!in2) w2=(LbmFloat)(sx)/(1.5*d2);
1175                                 if(!in1) w1=0.; if(!in2) w2=0.; // sharp falloff
1176                           uvec += v1*w1;
1177                           uvec += v2*w2;
1178                                 initVelocityCell(level, i,j,0, CFFluid, rho, rhomass, uvec );
1179                                 //errMsg("VORTT","init uvec"<<uvec);
1180                         }
1181
1182         }
1183 #endif // LBM_INCLUDE_TESTSOLVERS==1
1184
1185         //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
1186
1187         // prepare interface cells
1188         initFreeSurfaces();
1189         initStandingFluidGradient();
1190
1191         // perform first step to init initial mass
1192         mInitialMass = 0.0;
1193         int inmCellCnt = 0;
1194         FSGR_FORIJK1(mMaxRefine) {
1195                 if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
1196                         LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); 
1197                         FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
1198                         mInitialMass += fluidRho;
1199                         inmCellCnt ++;
1200                 } else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
1201                         mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
1202                         inmCellCnt ++;
1203                 }
1204         }
1205         mCurrentVolume = mCurrentMass = mInitialMass;
1206
1207         ParamVec cspv = mpParam->calculateCellSize();
1208         if(LBMDIM==2) cspv[2] = 1.0;
1209         inmCellCnt = 1;
1210         double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
1211         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
1212         mInitialMass = 0.0; // reset, and use actual value after first step
1213
1214         //mStartSymm = false;
1215 #if ELBEEM_PLUGIN!=1
1216         if((LBMDIM==2)&&(mSizex<200)) {
1217                 if(!checkSymmetry("init")) {
1218                         errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
1219                 } else {
1220                         errMsg("LbmFsgrSolver::initialize","Symmetric init!");
1221                 }
1222         }
1223 #endif // ELBEEM_PLUGIN!=1
1224         return true;
1225 }
1226
1227
1228 /*! prepare actual simulation start, setup viz etc */
1229 bool LbmFsgrSolver::initializeSolverPostinit() {
1230         // coarsen region
1231         myTime_t fsgrtstart = getTime(); 
1232         for(int lev=mMaxRefine-1; lev>=0; lev--) {
1233                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
1234                 adaptGrid(lev);
1235                 coarseRestrictFromFine(lev);
1236                 adaptGrid(lev);
1237                 coarseRestrictFromFine(lev);
1238         }
1239         markedClearList();
1240         myTime_t fsgrtend = getTime(); 
1241         if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
1242         mNumFsgrChanges = 0;
1243
1244         for(int l=0; l<cDirNum; l++) { 
1245                 LbmFloat area = 0.5 * 0.5 *0.5;
1246                 if(LBMDIM==2) area = 0.5 * 0.5;
1247
1248                 if(dfVecX[l]!=0) area *= 0.5;
1249                 if(dfVecY[l]!=0) area *= 0.5;
1250                 if(dfVecZ[l]!=0) area *= 0.5;
1251                 mFsgrCellArea[l] = area;
1252         } // l
1253
1254         // make sure both sets are ok
1255         // copy from other to curr
1256         for(int lev=0; lev<=mMaxRefine; lev++) {
1257         FSGR_FORIJK_BOUNDS(lev) {
1258                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1259         } } // first copy flags */
1260
1261
1262         // old mpPreviewSurface init
1263         //if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
1264         // make sure fill fracs are right for first surface generation
1265         stepMain();
1266
1267         // prepare once...
1268         mpIso->setParticles(mpParticles, mPartDropMassSub);
1269   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
1270         prepareVisualization();
1271         // copy again for stats counting
1272         for(int lev=0; lev<=mMaxRefine; lev++) {
1273         FSGR_FORIJK_BOUNDS(lev) {
1274                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1275         } } // first copy flags */
1276
1277
1278         // now really done...
1279   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
1280   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
1281         mInitDone = 1;
1282
1283         // init fluid control
1284         initCpdata();
1285
1286 #if LBM_INCLUDE_TESTSOLVERS==1
1287         initTestdata();
1288 #endif // ELBEEM_PLUGIN!=1
1289         // not inited? dont use...
1290         if(mCutoff<0) mCutoff=0;
1291
1292         initParticles();
1293         return true;
1294 }
1295
1296
1297
1298 // macros for mov obj init
1299 #if LBMDIM==2
1300
1301 #define POS2GRID_CHECK(vec,n) \
1302                                 monTotal++;\
1303                                 int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
1304                                 if(k!=0) continue; \
1305                                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
1306                                 if(i<=0) continue; \
1307                                 if(i>=mLevel[level].lSizex-1) continue; \
1308                                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
1309                                 if(j<=0) continue; \
1310                                 if(j>=mLevel[level].lSizey-1) continue;  \
1311
1312 #else // LBMDIM -> 3
1313 #define POS2GRID_CHECK(vec,n) \
1314                                 monTotal++;\
1315                                 const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
1316                                 if(i<=0) continue; \
1317                                 if(i>=mLevel[level].lSizex-1) continue; \
1318                                 const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
1319                                 if(j<=0) continue; \
1320                                 if(j>=mLevel[level].lSizey-1) continue; \
1321                                 const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
1322                                 if(k<=0) continue; \
1323                                 if(k>=mLevel[level].lSizez-1) continue; \
1324
1325 #endif // LBMDIM 
1326
1327 // calculate object velocity from vert arrays in objvel vec
1328 #define OBJVEL_CALC \
1329                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
1330                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
1331                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
1332                                 if(usqr>maxusqr) { \
1333                                         /* cutoff at maxVelVal */ \
1334                                         for(int jj=0; jj<3; jj++) { \
1335                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  \
1336                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
1337                                         } \
1338                                 } } \
1339                                 if(ntype&(CFBndFreeslip)) { \
1340                                         const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
1341                                         const LbmVec oldov=objvel; /*DEBUG*/ \
1342                                         objvel = vec2L((*pNormals)[n]) *dp; \
1343                                         /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
1344                                 } \
1345                                 else if(ntype&(CFBndPartslip)) { \
1346                                         const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
1347                                         const LbmVec oldov=objvel; /*DEBUG*/ \
1348                                         /* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
1349                                         const LbmFloat partv = mObjectPartslips[OId]; \
1350                                         /*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
1351                                         /* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
1352                                         objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
1353                                 }
1354
1355 #define TTT \
1356
1357
1358 /*****************************************************************************/
1359 //! init moving obstacles for next sim step sim 
1360 /*****************************************************************************/
1361 void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
1362         myTime_t monstart = getTime();
1363
1364         // movobj init
1365         const int level = mMaxRefine;
1366         const int workSet = mLevel[level].setCurr;
1367         const int otherSet = mLevel[level].setOther;
1368         LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
1369         // for debugging - check targetTime check during DEFAULT STREAM
1370         LbmFloat targetTime = mSimulationTime + mpParam->getTimestep();
1371         if(mLastSimTime == targetTime) {
1372                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
1373                 return;
1374         }
1375         //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
1376         //if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);
1377
1378         const LbmFloat maxVelVal = 0.1666;
1379         const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
1380
1381         LbmFloat rhomass = 0.0;
1382         CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
1383         CellFlagType ntype = CFInvalid;
1384         // WARNING - copied from geo init!
1385         int numobjs = (int)(mpGiObjects->size());
1386         ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1387         // 2d display as rectangles
1388         ntlVec3Gfx iniPos(0.0);
1389         if(LBMDIM==2) {
1390                 dvec[2] = 1.0; 
1391                 iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
1392         } else {
1393                 iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
1394         }
1395         
1396         if( (int)mObjectMassMovnd.size() < numobjs) {
1397                 for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
1398                         mObjectMassMovnd.push_back(0.);
1399                 }
1400         }
1401         
1402         // stats
1403         int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
1404         int nbored;
1405         for(int OId=0; OId<numobjs; OId++) {
1406                 ntlGeometryObject *obj = (*mpGiObjects)[OId];
1407                 bool skip = false;
1408                 if(obj->getGeoInitId() != mLbmInitId) skip=true;
1409                 if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
1410                 if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
1411                 if(skip) continue;
1412                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);
1413
1414                 if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || 
1415                                 ((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) {
1416
1417                         otype = ntype = CFInvalid;
1418                         switch(obj->getGeoInitType()) {
1419                                 /* case FGI_BNDPART: // old, use noslip for moving part/free objs
1420                                 case FGI_BNDFREE: 
1421                                         if(!staticInit) {
1422                                                 errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
1423                                                 otype = ntype = CFBnd|CFBndNoslip;
1424                                         } else {
1425                                                 if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
1426                                                 if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
1427                                         }
1428                                         break; 
1429                                         // off */
1430                                 case FGI_BNDPART: rhomass = BND_FILL;
1431                                         otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
1432                                         break;
1433                                 case FGI_BNDFREE: rhomass = BND_FILL;
1434                                         otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
1435                                         break;
1436                                         // off */
1437                                 case FGI_BNDNO:   rhomass = BND_FILL;
1438                                         otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
1439                                         break;
1440                                 case FGI_FLUID: 
1441                                         otype = ntype = CFFluid; 
1442                                         break;
1443                                 case FGI_MBNDINFLOW: 
1444                                         otype = ntype = CFMbndInflow; 
1445                                         break;
1446                                 case FGI_MBNDOUTFLOW: 
1447                                         otype = ntype = CFMbndOutflow; 
1448                                         break;
1449                         }
1450                         int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
1451                         int active =    ((obj->getGeoActive(targetTime)>0.)? 1:0);
1452                         //errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<"  s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
1453                         // skip inactive in/out flows
1454                         if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
1455                         if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;
1456
1457                         // copied from  recalculateObjectSpeeds
1458                         mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
1459                         debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
1460
1461                         //vector<ntlVec3Gfx> tNormals;
1462                         vector<ntlVec3Gfx> *pNormals = NULL;
1463                         mMOINormals.clear();
1464                         if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
1465
1466                         mMOIVertices.clear();
1467                         if(obj->getMeshAnimated()) { 
1468                                 // do two full update
1469                                 // TODO tNormals handling!?
1470                                 mMOIVerticesOld.clear();
1471                                 obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals,  mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
1472                                 monTrafo += mMOIVerticesOld.size();
1473                                 obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
1474                                 monTrafo += mMOIVertices.size();
1475                                 obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
1476                         } else {
1477                                 // only do transform update
1478                                 obj->getMovingPoints(mMOIVertices,pNormals); // mMOIVertices = mCachedMovPoints
1479                                 mMOIVerticesOld = mMOIVertices;
1480                                 // WARNING - assumes mSimulationTime is global!?
1481                                 obj->applyTransformation(targetTime, &mMOIVertices,pNormals, 0, mMOIVertices.size(), false );
1482                                 monTrafo += mMOIVertices.size();
1483
1484                                 // correct flags from last position, but extrapolate
1485                                 // velocity to next timestep
1486                                 obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
1487                                 monTrafo += mMOIVerticesOld.size();
1488                         }
1489
1490                         // object types
1491                         if(ntype&CFBnd){
1492
1493                                 // check if object is moving at all
1494                                 if(obj->getIsAnimated()) {
1495                                         ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
1496                                         // FIXME?
1497                                         if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
1498                                         // get old type - CHECK FIXME , timestep could have changed - cause trouble?
1499                                         ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
1500                                         if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
1501                                 }
1502                                 if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
1503                                 CellFlagType rflagnb[27];
1504                                 LbmFloat massCheck = 0.;
1505                                 int massReinits=0;                              
1506                                 bool fillCells = (mObjectMassMovnd[OId]<=-1.);
1507                                 LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
1508                                 
1509
1510                                 // first pass - set new obs. cells
1511                                 if(active) {
1512                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1513                                                 //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
1514                                                 POS2GRID_CHECK(mMOIVertices,n);
1515                                                 //{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
1516                                                 if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
1517                                                 monPoints++;
1518                                                 
1519                                                 // check mass
1520                                                 if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
1521                                                         FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
1522                                                         massReinits++;
1523                                                 }
1524                                                 else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
1525                                                         massCheck -= QCELL(level, i,j,k, workSet, dMass);
1526                                                         massReinits++;
1527                                                 }
1528
1529                                                 RFLAG(level, i,j,k, workSet) = ntype;
1530                                                 FORDF1 {
1531                                                         //CellFlagType flag = RFLAG_NB(level, i,j,k,workSet,l);
1532                                                         rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
1533                                                         if(rflagnb[l]&(CFFluid|CFInter)) {
1534                                                                 rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
1535                                                                 RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l]; 
1536                                                         }
1537                                                 }
1538                                                 LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
1539                                                 RAC(dstCell,0) = 0.0;
1540                                                 if(ntype&CFBndMoving) {
1541                                                         OBJVEL_CALC;
1542                                                         
1543                                                         // compute fluid acceleration
1544                                                         FORDF1 {
1545                                                                 if(rflagnb[l]&(CFFluid|CFInter)) { 
1546                                                                         const LbmFloat ux = dfDvecX[l]*objvel[0];
1547                                                                         const LbmFloat uy = dfDvecY[l]*objvel[1];
1548                                                                         const LbmFloat uz = dfDvecZ[l]*objvel[2];
1549
1550                                                                         LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // 
1551                                                                         if(ntype&(CFBndFreeslip|CFBndPartslip)) {
1552                                                                                 // missing, diag mass checks...
1553                                                                                 //if(l<=LBMDIM*2) factor *= 1.4142;
1554                                                                                 factor *= 2.0; // TODO, optimize
1555                                                                         } else {
1556                                                                                 factor *= 1.2; // TODO, optimize
1557                                                                         }
1558                                                                         factor *= impactCorrFactor; // use manual tweaking channel
1559                                                                         RAC(dstCell,l) = factor;
1560                                                                         massCheck += factor;
1561                                                                 } else {
1562                                                                         RAC(dstCell,l) = 0.;
1563                                                                 }
1564                                                         }
1565
1566 #if NEWDIRVELMOTEST==1
1567                                                         FORDF1 { RAC(dstCell,l) = 0.; }
1568                                                         RAC(dstCell,dMass)  = objvel[0];
1569                                                         RAC(dstCell,dFfrac) = objvel[1];
1570                                                         RAC(dstCell,dC)     = objvel[2];
1571 #endif // NEWDIRVELMOTEST==1
1572                                                 } else {
1573                                                         FORDF1 { RAC(dstCell,l) = 0.0; }
1574                                                 }
1575                                                 RAC(dstCell, dFlux) = targetTime;
1576                                                 //errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" dflt="<<RAC(dstCell, dFlux) );
1577                                                 monObsts++;
1578                                         } // points
1579                                 } // bnd, is active?
1580
1581                                 // second pass, remove old ones
1582                                 // warning - initEmptyCell et. al dont overwrite OId or persist flags...
1583                                 if(wasActive) {
1584                                         for(size_t n=0; n<mMOIVerticesOld.size(); n++) {
1585                                                 POS2GRID_CHECK(mMOIVerticesOld,n);
1586                                                 monPoints++;
1587                                                 if((RFLAG(level, i,j,k, workSet) == otype) &&
1588                                                          (QCELL(level, i,j,k, workSet, dFlux) != targetTime)) {
1589                                                         // from mainloop
1590                                                         nbored = 0;
1591                                                         // TODO: optimize for OPT3D==0
1592                                                         FORDF1 {
1593                                                                 //rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
1594                                                                 rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance
1595                                                                 nbored |= rflagnb[l];
1596                                                         } 
1597                                                         CellFlagType settype = CFInvalid;
1598                                                         if(nbored&CFFluid) {
1599                                                                 settype = CFInter|CFNoInterpolSrc; 
1600                                                                 rhomass = 1.5;
1601                                                                 if(!fillCells) rhomass = 0.;
1602
1603                                                                 OBJVEL_CALC;
1604                                                                 if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
1605
1606                                                                 // new interpolate values
1607                                                                 LbmFloat avgrho = 0.0;
1608                                                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1609                                                                 interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
1610                                                                 initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
1611                                                                 //errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
1612                                                                 massCheck += rhomass;
1613                                                         } else {
1614                                                                 settype = CFEmpty; rhomass = 0.0;
1615                                                                 initEmptyCell(level, i,j,k, settype, 1., rhomass );
1616                                                         }
1617                                                         monFluids++;
1618                                                         massReinits++;
1619                                                 } // flag & simtime
1620                                         }
1621                                 }  // wasactive
1622
1623                                 // only compute mass transfer when init is done
1624                                 if(mInitDone) {
1625                                         errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
1626                                                 " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass ); 
1627                                         mObjectMassMovnd[OId] += massCheck;
1628                                 }
1629                         } // bnd, active
1630
1631                         else if(ntype&CFFluid){
1632                                 // second static init pass
1633                                 if(staticInit) {
1634                                         //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" verts"<<mMOIVertices.size() ,9);
1635                                         CellFlagType setflflag = CFFluid; //|(OId<<24);
1636                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1637                                                 POS2GRID_CHECK(mMOIVertices,n);
1638                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1639                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1640                                                         //changeFlag(level, i,j,k, workSet, setflflag);
1641                                                         initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
1642                                                 } 
1643                                                 //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
1644                                         }
1645                                 } // second static inflow pass
1646                         } // fluid
1647
1648                         else if(ntype&CFMbndInflow){
1649                                 // inflow pass - add new fluid cells
1650                                 // this is slightly different for standing inflows,
1651                                 // as the fluid is forced to the desired velocity inside as 
1652                                 // well...
1653                                 const LbmFloat iniRho = 1.0;
1654                                 const LbmVec vel(mObjectSpeeds[OId]);
1655                                 const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5;
1656                                 USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
1657                                 //errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
1658                                 
1659                                 for(size_t n=0; n<mMOIVertices.size(); n++) {
1660                                         POS2GRID_CHECK(mMOIVertices,n);
1661                                         // TODO - also reinit interface cells !?
1662                                         if(RFLAG(level, i,j,k, workSet)!=CFEmpty) { 
1663                                                 // test prevent particle gen for inflow inter cells
1664                                                 if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag
1665                                                 continue; }
1666                                         monFluids++;
1667
1668                                         // TODO add OPT3D treatment
1669                                         LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
1670                                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
1671                                         RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
1672                                         RAC(tcel, dFlux) = FLUX_INIT;
1673                                         CellFlagType setFlag = CFInter;
1674                                         changeFlag(level, i,j,k, workSet, setFlag);
1675                                         mInitialMass += iniRho;
1676                                 }
1677                                 // second static init pass
1678                                 if(staticInit) {
1679                                         CellFlagType set2Flag = CFMbndInflow|(OId<<24);
1680                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1681                                                 POS2GRID_CHECK(mMOIVertices,n);
1682                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1683                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1684                                                         forceChangeFlag(level, i,j,k, workSet, set2Flag);
1685                                                 } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
1686                                                         forceChangeFlag(level, i,j,k, workSet, 
1687                                                                         (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
1688                                                 }
1689                                         }
1690                                 } // second static inflow pass
1691
1692                         } // inflow
1693
1694                         else if(ntype&CFMbndOutflow){
1695                                 const LbmFloat iniRho = 0.0;
1696                                 for(size_t n=0; n<mMOIVertices.size(); n++) {
1697                                         POS2GRID_CHECK(mMOIVertices,n);
1698                                         // FIXME check fluid/inter cells for non-static!?
1699                                         if(!(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter))) {
1700                                                 if((staticInit)&&(RFLAG(level, i,j,k, workSet)==CFEmpty)) {
1701                                                         forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
1702                                                 continue;
1703                                         }
1704                                         monFluids++;
1705                                         // interface cell - might be double empty? FIXME check?
1706
1707                                         // remove fluid cells, shouldnt be here anyway
1708                                         LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
1709                                         LbmFloat fluidRho = RAC(tcel,0); FORDF1 { fluidRho += RAC(tcel,l); }
1710                                         mInitialMass -= fluidRho;
1711                                         RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
1712                                         RAC(tcel, dFlux) = FLUX_INIT;
1713                                         CellFlagType setFlag = CFInter;
1714                                         changeFlag(level, i,j,k, workSet, setFlag);
1715
1716                                         // same as ifemptied for if below
1717                                         LbmPoint emptyp;
1718                                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1719                                         mListEmpty.push_back( emptyp );
1720                                         //calcCellsEmptied++;
1721                                 } // points
1722                                 // second static init pass
1723                                 if(staticInit) {
1724                                         CellFlagType set2Flag = CFMbndOutflow|(OId<<24);
1725                                         for(size_t n=0; n<mMOIVertices.size(); n++) {
1726                                                 POS2GRID_CHECK(mMOIVertices,n);
1727                                                 if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
1728                                                 if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
1729                                                         forceChangeFlag(level, i,j,k, workSet, set2Flag);
1730                                                 } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
1731                                                         forceChangeFlag(level, i,j,k, workSet, 
1732                                                                         (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
1733                                                 }
1734                                         }
1735                                 } // second static outflow pass
1736                         } // outflow
1737
1738                 } // allbound check
1739         } // OId
1740
1741
1742         /* { // DEBUG check
1743                 int workSet = mLevel[level].setCurr;
1744                 FSGR_FORIJK1(level) {
1745                         if( (RFLAG(level,i,j,k, workSet) & CFBndMoving) ) {
1746                                 if(QCELL(level, i,j,k, workSet, dFlux)!=targetTime) {
1747                                         errMsg("lastt"," old val!? at "<<PRINT_IJK<<" t="<<QCELL(level, i,j,k, workSet, dFlux)<<" target="<<targetTime);
1748                                 }
1749                         }
1750                 }
1751         } // DEBUG */
1752         
1753 #undef POS2GRID_CHECK
1754         myTime_t monend = getTime();
1755         if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
1756         mLastSimTime = targetTime;
1757 }
1758
1759
1760 // geoinit position
1761
1762 #define GETPOS(i,j,k)  \
1763         ntlVec3Gfx ggpos = \
1764                 ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
1765                                   iniPos[1]+ dvec[1]*(gfxReal)(j), \
1766                                   iniPos[2]+ dvec[2]*(gfxReal)(k) ); \
1767   if((LBMDIM==2)&&(mInit2dYZ)) { SWAPYZ(ggpos); }
1768
1769 /*****************************************************************************/
1770 /*! perform geometry init (if switched on) */
1771 /*****************************************************************************/
1772 extern int globGeoInitDebug; //solver_interface
1773 bool LbmFsgrSolver::initGeometryFlags() {
1774         int level = mMaxRefine;
1775         myTime_t geotimestart = getTime(); 
1776         ntlGeometryObject *pObj;
1777         ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1778         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
1779         // WARNING - copied to movobj init!
1780
1781         /* init object velocities, this has always to be called for init */
1782         initGeoTree();
1783         if(mAllfluid) { 
1784                 freeGeoTree();
1785                 return true; }
1786
1787         // make sure moving obstacles are inited correctly
1788         // preinit moving obj points
1789         int numobjs = (int)(mpGiObjects->size());
1790         for(int o=0; o<numobjs; o++) {
1791                 ntlGeometryObject *obj = (*mpGiObjects)[o];
1792                 //debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
1793                 if(
1794                                 ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
1795                                 (obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
1796                         if(!obj->getMeshAnimated()) { 
1797                                 debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
1798                                 obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
1799                         }
1800                 }
1801         }
1802
1803         // max speed init
1804         ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
1805         ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
1806         mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
1807         LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
1808         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
1809         if(mpParam->getSimulationMaxSpeed() > allowMax) {
1810                 // similar to adaptTimestep();
1811                 LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
1812                 LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
1813                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
1814                 mpParam->setDesiredTimestep( newdt );
1815                 mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
1816                 maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
1817                                       mSimulationTime, mpParam->getTimestep() )) ));
1818                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
1819         }
1820         recalculateObjectSpeeds();
1821         // */
1822
1823         // init obstacles for first time step (requires obj speeds)
1824         initMovingObstacles(true);
1825
1826         /* set interface cells */
1827         ntlVec3Gfx pos,iniPos; // position of current cell
1828         LbmFloat rhomass = 0.0;
1829         CellFlagType ntype = CFInvalid;
1830         int savedNodes = 0;
1831         int OId = -1;
1832         gfxReal distance;
1833
1834         // 2d display as rectangles
1835         if(LBMDIM==2) {
1836                 dvec[2] = 0.0; 
1837                 iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
1838                 //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
1839         } else {
1840                 iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
1841         }
1842
1843
1844         // first init boundary conditions
1845         // invalid cells are set to empty afterwards
1846         // TODO use floop macros!?
1847         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1848                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1849                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1850                                 ntype = CFInvalid;
1851                                 
1852                                 GETPOS(i,j,k);
1853                                 const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
1854                                 if(inside) {
1855                                         pObj = (*mpGiObjects)[OId];
1856                                         switch( pObj->getGeoInitType() ){
1857                                         case FGI_MBNDINFLOW:  
1858                                           if(! pObj->getIsAnimated() ) {
1859                                                         rhomass = 1.0;
1860                                                         ntype = CFFluid | CFMbndInflow;
1861                                                 } else {
1862                                                         ntype = CFInvalid;
1863                                                 }
1864                                                 break;
1865                                         case FGI_MBNDOUTFLOW: 
1866                                           if(! pObj->getIsAnimated() ) {
1867                                                         rhomass = 0.0;
1868                                                         ntype = CFEmpty|CFMbndOutflow; 
1869                                                 } else {
1870                                                         ntype = CFInvalid;
1871                                                 }
1872                                                 break;
1873                                         case FGI_BNDNO: 
1874                                                 rhomass = BND_FILL;
1875                                                 ntype = CFBnd|CFBndNoslip; 
1876                                                 break;
1877                                         case FGI_BNDPART: 
1878                                                 rhomass = BND_FILL;
1879                                                 ntype = CFBnd|CFBndPartslip; break;
1880                                         case FGI_BNDFREE: 
1881                                                 rhomass = BND_FILL;
1882                                                 ntype = CFBnd|CFBndFreeslip; break;
1883                                         default: // warn here?
1884                                                 rhomass = BND_FILL;
1885                                                 ntype = CFBnd|CFBndNoslip; break;
1886                                         }
1887                                 }
1888                                 if(ntype != CFInvalid) {
1889                                         // initDefaultCell
1890                                         if((ntype & CFMbndInflow) || (ntype & CFMbndOutflow) ) { }
1891                                         ntype |= (OId<<24); // NEWTEST2
1892                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1893                                 }
1894
1895                                 // walk along x until hit for following inits
1896                                 if(distance<=-1.0) { distance = 100.0; } // FIXME dangerous
1897                                 if(distance>=0.0) {
1898                                         gfxReal dcnt=dvec[0];
1899                                         while(( dcnt< distance-dvec[0] )&&(i+1<mLevel[level].lSizex-1)) {
1900                                                 dcnt += dvec[0]; i++;
1901                                                 savedNodes++;
1902                                                 if(ntype != CFInvalid) {
1903                                                         // rho,mass,OId are still inited from above
1904                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1905                                                 }
1906                                         }
1907                                 } 
1908                                 // */
1909
1910                         } 
1911                 } 
1912         } // zmax
1913         // */
1914
1915         // now init fluid layer
1916         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1917                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1918                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1919                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
1920                                 ntype = CFInvalid;
1921                                 int inits = 0;
1922                                 GETPOS(i,j,k);
1923                                 const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
1924                                 if(inside) {
1925                                         ntype = CFFluid;
1926                                 }
1927                                 if(ntype != CFInvalid) {
1928                                         // initDefaultCell
1929                                         rhomass = 1.0;
1930                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1931                                         inits++;
1932                                 }
1933
1934                                 // walk along x until hit for following inits
1935                                 if(distance<=-1.0) { distance = 100.0; }
1936                                 if(distance>=0.0) {
1937                                         gfxReal dcnt=dvec[0];
1938                                         while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
1939                                                 dcnt += dvec[0]; i++;
1940                                                 savedNodes++;
1941                                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
1942                                                 if(ntype != CFInvalid) {
1943                                                         // rhomass are still inited from above
1944                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
1945                                                         inits++;
1946                                                 }
1947                                         }
1948                                 } // distance>0
1949                                 
1950                         } 
1951                 } 
1952         } // zmax
1953
1954         // reset invalid to empty again
1955         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
1956                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
1957                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
1958                                 if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) {
1959                                         RFLAG(level, i,j,k, mLevel[level].setOther) = 
1960                                         RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty;
1961                                 }
1962                         }
1963                 }
1964         }
1965
1966         freeGeoTree();
1967         myTime_t geotimeend = getTime(); 
1968         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 ); 
1969         //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
1970         return true;
1971 }
1972
1973 #undef GETPOS
1974
1975 /*****************************************************************************/
1976 /* init part for all freesurface testcases */
1977 void LbmFsgrSolver::initFreeSurfaces() {
1978         double interfaceFill = 0.45;   // filling level of interface cells
1979         //interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!
1980
1981         // set interface cells 
1982         FSGR_FORIJK1(mMaxRefine) {
1983                 if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
1984                         FORDF1 {
1985                                 int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
1986                                 if( ( RFLAG(mMaxRefine, ni, nj, nk,  mLevel[mMaxRefine].setCurr)& CFEmpty ) ) {
1987                                         LbmFloat arho=0., aux=0., auy=0., auz=0.;
1988                                         interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
1989                                         //errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
1990                                         // unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
1991                                         initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
1992                                 }
1993                         }
1994                 }
1995         }
1996
1997         // remove invalid interface cells 
1998         FSGR_FORIJK1(mMaxRefine) {
1999                 // remove invalid interface cells 
2000                 if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
2001                         int delit = 0;
2002                         int NBs = 0; // neighbor flags or'ed 
2003                         int noEmptyNB = 1;
2004
2005                         FORDF1 {
2006                                 if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
2007                                         noEmptyNB = 0;
2008                                 }
2009                                 NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
2010                         }
2011                         // remove cells with no fluid or interface neighbors
2012                         if((NBs & CFFluid)==0) { delit = 1; }
2013                         if((NBs & CFInter)==0) { delit = 1; }
2014                         // remove cells with no empty neighbors
2015                         if(noEmptyNB) { delit = 2; }
2016
2017                         // now we can remove the cell 
2018                         if(delit==1) {
2019                                 initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
2020                         }
2021                         if(delit==2) {
2022                                 //initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
2023                                 LbmFloat arho=0., aux=0., auy=0., auz=0.;
2024                                 interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
2025                                 initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
2026                         }
2027                 } // interface 
2028         } // */
2029
2030         // another brute force init, make sure the fill values are right...
2031         // and make sure both sets are equal
2032         for(int lev=0; lev<=mMaxRefine; lev++) {
2033         FSGR_FORIJK_BOUNDS(lev) {
2034                 if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { 
2035                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
2036                         continue;
2037                 }
2038                 if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { 
2039                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
2040                         continue;
2041                 }
2042         } }
2043
2044         // ----------------------------------------------------------------------
2045         // smoother surface...
2046         if(mInitSurfaceSmoothing>0) {
2047                 debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing timestep ",10);
2048 #if COMPRESSGRIDS==1
2049                 //errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
2050 #endif // COMPRESSGRIDS==0
2051         }
2052         for(int s=0; s<mInitSurfaceSmoothing; s++) {
2053                 //SGR_FORIJK1(mMaxRefine) {
2054
2055                 int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
2056                 int lev = mMaxRefine;
2057 #if COMPRESSGRIDS==0
2058                 for(int k=kstart;k<kend;++k) {
2059 #else // COMPRESSGRIDS==0
2060                 int kdir = 1; // COMPRT ON
2061                 if(mLevel[lev].setCurr==1) {
2062                         kdir = -1;
2063                         int temp = kend;
2064                         kend = kstart-1;
2065                         kstart = temp-1;
2066                 } // COMPRT
2067                 for(int k=kstart;k!=kend;k+=kdir) {
2068 #endif // COMPRESSGRIDS==0
2069                 for(int j=1;j<mLevel[lev].lSizey-1;++j) {
2070                 for(int i=1;i<mLevel[lev].lSizex-1;++i) {
2071                         if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
2072                                 LbmFloat mass = 0.0;
2073                                 //LbmFloat nbdiv;
2074                                 //FORDF0 {
2075                                 for(int l=0;(l<cDirNum); l++) { 
2076                                         int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
2077                                         if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
2078                                                 mass += 1.0;
2079                                         }
2080                                         if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFInter ){
2081                                                 mass += QCELL(lev, ni,nj,nk, mLevel[lev].setCurr, dMass);
2082                                         }
2083                                         //nbdiv+=1.0;
2084                                 }
2085
2086                                 //errMsg(" I ", PRINT_IJK<<" m"<<mass );
2087                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
2088                                 QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
2089                         }
2090                 }}}
2091
2092                 mLevel[lev].setOther = mLevel[lev].setCurr;
2093                 mLevel[lev].setCurr ^= 1;
2094         }
2095         // copy back...?
2096 }
2097
2098 /*****************************************************************************/
2099 /* init part for all freesurface testcases */
2100 void LbmFsgrSolver::initStandingFluidGradient() {
2101         // ----------------------------------------------------------------------
2102         // standing fluid preinit
2103         const int debugStandingPreinit = 0;
2104         int haveStandingFluid = 0;
2105
2106 #define STANDFLAGCHECK(iindex) \
2107                                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
2108                                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
2109                                         if((iindex)>1) { haveStandingFluid=(iindex); } \
2110                                         j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
2111                                         continue; \
2112                                 } 
2113         int gravIndex[3] = {0,0,0};
2114         int gravDir[3] = {1,1,1};
2115         int maxGravComp = 1; // by default y
2116         int gravComp1 = 0; // by default y
2117         int gravComp2 = 2; // by default y
2118         if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
2119         if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
2120
2121         int gravIMin[3] = { 0 , 0 , 0 };
2122         int gravIMax[3] = {
2123                 mLevel[mMaxRefine].lSizex + 0,
2124                 mLevel[mMaxRefine].lSizey + 0,
2125                 mLevel[mMaxRefine].lSizez + 0 };
2126         if(LBMDIM==2) gravIMax[2] = 1;
2127
2128         //int gravDir = 1;
2129         if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
2130                 // swap directions
2131                 int i=maxGravComp;
2132                 int tmp = gravIMin[i];
2133                 gravIMin[i] = gravIMax[i] - 1;
2134                 gravIMax[i] = tmp - 1;
2135                 gravDir[i] = -1;
2136         }
2137 #define PRINTGDIRS \
2138         errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
2139         errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
2140         errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); 
2141         // _PRINTGDIRS;
2142
2143         bool gravAbort = false;
2144 #define GRAVLOOP \
2145         gravAbort=false; \
2146         for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
2147                 for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
2148                         for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0]) 
2149
2150         GRAVLOOP {
2151                 int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2152                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || 
2153                     ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) || 
2154                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) {  
2155                         int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
2156                         if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
2157                         if(fluidHeight>1) {
2158                                 haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; 
2159                                 gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
2160                         }
2161                         gravAbort = true; continue; 
2162                 } 
2163         } // GRAVLOOP
2164         // _PRINTGDIRS;
2165
2166         LbmFloat fluidHeight;
2167         fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
2168 #if LBM_INCLUDE_TESTSOLVERS==1
2169         if(mUseTestdata) mpTest->mFluidHeight = (int)fluidHeight;
2170 #endif // ELBEEM_PLUGIN!=1
2171         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])<<
2172                         " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
2173                                 
2174         if(mDisableStandingFluidInit) {
2175                 debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
2176                 haveStandingFluid=0;
2177         }
2178
2179         // copy flags and init , as no flags will be changed during grav init
2180         // also important for corasening later on
2181         const int lev = mMaxRefine;
2182         CellFlagType nbflag[LBM_DFNUM], nbored; 
2183         for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
2184                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2185                         for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2186                                 if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
2187                                         nbored = 0;
2188                                         FORDF1 {
2189                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2190                                                 nbored |= nbflag[l];
2191                                         } 
2192                                         if(nbored&CFBnd) {
2193                                                 RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
2194                                         } else {
2195                                                 RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
2196                                         }
2197                                 }
2198                                 RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
2199         } } }
2200
2201         if(haveStandingFluid) {
2202                 int rhoworkSet = mLevel[lev].setCurr;
2203                 myTime_t timestart = getTime(); // FIXME use user time here?
2204
2205                 GRAVLOOP {
2206                         int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2207                         //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
2208                         if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
2209                                         ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ 
2210                                 //gravAbort = true; 
2211                                 continue;
2212                         }
2213
2214                         LbmFloat rho = 1.0;
2215                         // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
2216                         rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * 
2217                                 (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
2218                         if(debugStandingPreinit) 
2219                                 if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { 
2220                                         errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); 
2221                                 }
2222
2223                         if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
2224                                         (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
2225                                 FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
2226                                 QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
2227                         }
2228
2229                 } // GRAVLOOP
2230
2231                 debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<<
2232                         (1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8);
2233                 
2234                 int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
2235                 preinitSteps = (haveStandingFluid>>2); // not much use...?
2236                 //preinitSteps = 0;
2237                 debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
2238                 for(int s=0; s<preinitSteps; s++) {
2239                         // in solver main cpp
2240                         standingFluidPreinit();
2241                 }
2242
2243                 myTime_t timeend = getTime();
2244                 debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
2245 #undef  NBFLAG
2246         }
2247 }
2248
2249
2250
2251 bool LbmFsgrSolver::checkSymmetry(string idstring)
2252 {
2253         bool erro = false;
2254         bool symm = true;
2255         int msgs = 0;
2256         const int maxMsgs = 10;
2257         const bool markCells = false;
2258
2259         //for(int lev=0; lev<=mMaxRefine; lev++) {
2260         { int lev = mMaxRefine;
2261
2262         // no point if not symm.
2263         if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
2264                 // ok
2265         } else {
2266                 return false;
2267         }
2268
2269         for(int s=0; s<2; s++) {
2270         FSGR_FORIJK1(lev) {
2271                 if(i<(mLevel[lev].lSizex/2)) {
2272                         int inb = (mLevel[lev].lSizey-1-i); 
2273
2274                         if(lev==mMaxRefine) inb -= 1;           // FSGR_SYMM_T
2275
2276                         if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
2277                                 if(LBMDIM==2) {
2278                                         if(msgs<maxMsgs) { msgs++;
2279                                                 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) );
2280                                         }
2281                                 }
2282                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2283                                 symm = false;
2284                         }
2285                         if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
2286                                 if(LBMDIM==2) {
2287                                         if(msgs<maxMsgs) { msgs++;
2288                                                 //debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
2289                                                 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) );
2290                                         }
2291                                 }
2292                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2293                                 symm = false;
2294                         }
2295
2296                         LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
2297                         FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
2298                         LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
2299                         FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
2300                         if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
2301                                 if(LBMDIM==2) {
2302                                         if(msgs<maxMsgs) { msgs++;
2303                                                 //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
2304                                                 errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
2305                                         }
2306                                 }
2307                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2308                                 symm = false;
2309                         }
2310                 }
2311         } }
2312         } // lev
2313         LbmFloat maxdiv =0.0;
2314         if(erro) {
2315                 errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
2316                 //if(LBMDIM==2) mPanic = true; 
2317                 //return false;
2318         } else {
2319                 errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
2320         }
2321         // all ok...
2322         return symm;
2323 }// */
2324
2325
2326 void 
2327 LbmFsgrSolver::interpolateCellValues(
2328                 int level,int ei,int ej,int ek,int workSet, 
2329                 LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz) 
2330 {
2331         LbmFloat avgrho = 0.0;
2332         LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
2333         LbmFloat cellcnt = 0.0;
2334
2335         for(int nbl=1; nbl< cDfNum ; ++nbl) {
2336                 if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
2337                 if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){
2338                                 //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
2339                                  //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) { 
2340                         cellcnt += 1.0;
2341                         for(int rl=0; rl< cDfNum ; ++rl) { 
2342                                 LbmFloat nbdf =  QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
2343                                 avgux  += (dfDvecX[rl]*nbdf); 
2344                                 avguy  += (dfDvecY[rl]*nbdf);  
2345                                 avguz  += (dfDvecZ[rl]*nbdf);  
2346                                 avgrho += nbdf;
2347                         }
2348                 }
2349         }
2350
2351         if(cellcnt<=0.0) {
2352                 // no nbs? just use eq.
2353                 avgrho = 1.0;
2354                 avgux = avguy = avguz = 0.0;
2355                 //TTT mNumProblems++;
2356 #if ELBEEM_PLUGIN!=1
2357                 //mPanic=1; 
2358                 // this can happen for worst case moving obj scenarios...
2359                 errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
2360 #endif // ELBEEM_PLUGIN
2361         } else {
2362                 // init speed
2363                 avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
2364                 avgrho /= cellcnt;
2365         }
2366
2367         retrho = avgrho;
2368         retux = avgux;
2369         retuy = avguy;
2370         retuz = avguz;
2371 } // interpolateCellValues
2372
2373
2374 /******************************************************************************
2375  * instantiation
2376  *****************************************************************************/
2377
2378 //! lbm factory functions
2379 LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }
2380
2381