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