← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/blender/gutsy-security

« back to all changes in this revision

Viewing changes to intern/elbeem/intern/solver_util.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Lukas Fittl
  • Date: 2006-09-20 01:57:27 UTC
  • mfrom: (1.2.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20060920015727-gmoqlxwstx9wwqs3
Tags: 2.42a-1ubuntu1
* Merge from Debian unstable (Closes: Malone #55903). Remaining changes:
  - debian/genpot: Add python scripts from Lee June <blender@eyou.com> to
    generate a reasonable PO template from the sources. Since gettext is used
    in a highly nonstandard way, xgettext does not work for this job.
  - debian/rules: Call the scripts, generate po/blender.pot, and clean it up
    in the clean target.
  - Add a proper header to the generated PO template.
* debian/control: Build depend on libavformat-dev >= 3:0.cvs20060823-3.1,
  otherwise this package will FTBFS

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/elbeem/intern/solver_util.cpp
7
7
 *
8
8
 *****************************************************************************/
9
9
 
10
 
#if ((!defined(__APPLE_CC__)) && (!defined(__INTEL_COMPILER))) || defined(LBM_FORCEINCLUDE)
11
10
#include "solver_class.h"
 
11
#include "solver_relax.h"
 
12
#include "particletracer.h"
12
13
 
13
14
 
14
15
/******************************************************************************
16
17
 *****************************************************************************/
17
18
 
18
19
//! for raytracing
19
 
template<class D>
20
 
void LbmFsgrSolver<D>::prepareVisualization( void ) {
 
20
void LbmFsgrSolver::prepareVisualization( void ) {
21
21
        int lev = mMaxRefine;
22
22
        int workSet = mLevel[lev].setCurr;
23
23
 
31
31
        for(int k= 0; k< 5; ++k) 
32
32
   for(int j=0;j<mLevel[lev].lSizey-0;j++) 
33
33
    for(int i=0;i<mLevel[lev].lSizex-0;i++) {
34
 
                *D::mpIso->lbmGetData(i,j,ZKOFF)=0.0;
 
34
                *this->mpIso->lbmGetData(i,j,ZKOFF)=0.0;
35
35
        }
36
36
#else // LBMDIM==2
37
37
        // 3d, use normal bounds
41
41
        for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
42
42
   for(int j=0;j<mLevel[lev].lSizey-0;j++) 
43
43
    for(int i=0;i<mLevel[lev].lSizex-0;i++) {
44
 
                *D::mpIso->lbmGetData(i,j,ZKOFF)=0.0;
 
44
                *this->mpIso->lbmGetData(i,j,ZKOFF)=0.0;
45
45
        }
46
46
#endif // LBMDIM==2
47
47
 
48
48
 
49
 
        
 
49
        LbmFloat minval = this->mIsoValue*1.1; // / mIsoWeight[13]; 
50
50
        // add up...
51
51
        float val = 0.0;
52
52
        for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 
53
53
   for(int j=1;j<mLevel[lev].lSizey-1;j++) 
54
54
    for(int i=1;i<mLevel[lev].lSizex-1;i++) {
55
 
 
56
 
                //continue; // OFF DEBUG
57
 
                if(RFLAG(lev, i,j,k,workSet)&(CFBnd|CFEmpty)) {
58
 
                        continue;
59
 
                } else
60
 
                if( (RFLAG(lev, i,j,k,workSet)&CFInter) && (!(RFLAG(lev, i,j,k,workSet)&CFNoNbEmpty)) ){
61
 
                        // no empty nb interface cells are treated as full
62
 
                        val =  (QCELL(lev, i,j,k,workSet, dFfrac)); 
63
 
                        /* // flicker-test-fix: no real difference
64
 
                        if( (!(RFLAG(lev, i,j,k,workSet)&CFNoBndFluid)) && 
65
 
                                                (RFLAG(lev, i,j,k,workSet)&CFNoNbFluid)   &&
66
 
                                                (val<D::mIsoValue) ){ 
67
 
                                        val = D::mIsoValue*1.1; }
68
 
                        // */
69
 
                } else {
70
 
                        // fluid?
71
 
                        val = 1.0; ///27.0;
72
 
                } // */
73
 
 
74
 
                *D::mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
75
 
                *D::mpIso->lbmGetData( i   , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 
76
 
                *D::mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); 
77
 
                                                                                
78
 
                *D::mpIso->lbmGetData( i-1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); 
79
 
                *D::mpIso->lbmGetData( i   , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); 
80
 
                *D::mpIso->lbmGetData( i+1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); 
81
 
                                                                                
82
 
                *D::mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); 
83
 
                *D::mpIso->lbmGetData( i   , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); 
84
 
                *D::mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); 
85
 
                                                                                
86
 
                                                                                
87
 
                *D::mpIso->lbmGetData( i-1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[9] ); 
88
 
                *D::mpIso->lbmGetData( i   , j-1  ,ZKOFF  ) += ( val * mIsoWeight[10] ); 
89
 
                *D::mpIso->lbmGetData( i+1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[11] ); 
90
 
                                                                                                                                        
91
 
                *D::mpIso->lbmGetData( i-1 , j    ,ZKOFF  ) += ( val * mIsoWeight[12] ); 
92
 
                *D::mpIso->lbmGetData( i   , j    ,ZKOFF  ) += ( val * mIsoWeight[13] ); 
93
 
                *D::mpIso->lbmGetData( i+1 , j    ,ZKOFF  ) += ( val * mIsoWeight[14] ); 
94
 
                                                                                                                                        
95
 
                *D::mpIso->lbmGetData( i-1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[15] ); 
96
 
                *D::mpIso->lbmGetData( i   , j+1  ,ZKOFF  ) += ( val * mIsoWeight[16] ); 
97
 
                *D::mpIso->lbmGetData( i+1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[17] ); 
98
 
                                                                                
99
 
                                                                                
100
 
                *D::mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); 
101
 
                *D::mpIso->lbmGetData( i   , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); 
102
 
                *D::mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); 
103
 
                                                                                                                                 
104
 
                *D::mpIso->lbmGetData( i-1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); 
105
 
                *D::mpIso->lbmGetData( i   , j   ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); 
106
 
                *D::mpIso->lbmGetData( i+1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); 
107
 
                                                                                                                                 
108
 
                *D::mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); 
109
 
                *D::mpIso->lbmGetData( i   , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); 
110
 
                *D::mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
 
55
                        const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
 
56
                        //if(cflag&(CFBnd|CFEmpty)) {
 
57
                        if(cflag&(CFBnd)) {
 
58
                                continue;
 
59
 
 
60
                        } else if( (cflag&CFEmpty) ) {
 
61
                                //continue; // OFF DEBUG
 
62
                                int noslipbnd = 0;
 
63
                                int intercnt = 0;
 
64
                                CellFlagType nbored;
 
65
                                FORDF1 { 
 
66
                                        const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
 
67
                                        if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; }
 
68
                                        if(nbflag&CFInter){ intercnt++; }
 
69
                                        nbored |= nbflag;
 
70
                                }
 
71
                                //? val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
 
72
                                if((noslipbnd)&&(intercnt>6)) {
 
73
                                        //if(val<minval) val = minval; 
 
74
                                        //*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
 
75
                                        *this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
 
76
                                } else if((noslipbnd)&&(intercnt>0)) {
 
77
                                        *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.95;
 
78
                                } else {
 
79
                                }
 
80
                                continue;
 
81
 
 
82
                        //} else if( (cflag&CFInter) ) {
 
83
 
 
84
                        } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
 
85
                                // optimized fluid
 
86
                                val = 1.;
 
87
 
 
88
                        } else if( (cflag&(CFInter|CFFluid)) ) {
 
89
                                int noslipbnd = 0;
 
90
                                FORDF1 { 
 
91
                                        const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
 
92
                                        if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; l=100; } 
 
93
                                }
 
94
                                // no empty nb interface cells are treated as full
 
95
                                if(cflag&(CFNoNbEmpty|CFFluid)) {
 
96
                                        val=1.0;
 
97
                                }
 
98
                                val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
 
99
                                
 
100
                                if(noslipbnd) {
 
101
                                        //errMsg("NEWVAL", PRINT_IJK<<" val"<<val <<" f"<< convertCellFlagType2String(cflag)<<" "<<noslipbnd); //" nbfl"<<convertCellFlagType2String(nbored) );
 
102
                                        if(val<minval) val = minval; 
 
103
                                        *this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
 
104
                                }
 
105
                                // */
 
106
 
 
107
                        } else { // unused?
 
108
                                continue;
 
109
                                // old fluid val = 1.0; 
 
110
                        } // */
 
111
 
 
112
                        *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
 
113
                        *this->mpIso->lbmGetData( i   , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 
 
114
                        *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); 
 
115
 
 
116
                        *this->mpIso->lbmGetData( i-1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); 
 
117
                        *this->mpIso->lbmGetData( i   , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); 
 
118
                        *this->mpIso->lbmGetData( i+1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); 
 
119
 
 
120
                        *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); 
 
121
                        *this->mpIso->lbmGetData( i   , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); 
 
122
                        *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); 
 
123
 
 
124
 
 
125
                        *this->mpIso->lbmGetData( i-1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[9] ); 
 
126
                        *this->mpIso->lbmGetData( i   , j-1  ,ZKOFF  ) += ( val * mIsoWeight[10] ); 
 
127
                        *this->mpIso->lbmGetData( i+1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[11] ); 
 
128
 
 
129
                        *this->mpIso->lbmGetData( i-1 , j    ,ZKOFF  ) += ( val * mIsoWeight[12] ); 
 
130
                        *this->mpIso->lbmGetData( i   , j    ,ZKOFF  ) += ( val * mIsoWeight[13] ); 
 
131
                        *this->mpIso->lbmGetData( i+1 , j    ,ZKOFF  ) += ( val * mIsoWeight[14] ); 
 
132
 
 
133
                        *this->mpIso->lbmGetData( i-1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[15] ); 
 
134
                        *this->mpIso->lbmGetData( i   , j+1  ,ZKOFF  ) += ( val * mIsoWeight[16] ); 
 
135
                        *this->mpIso->lbmGetData( i+1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[17] ); 
 
136
 
 
137
 
 
138
                        *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); 
 
139
                        *this->mpIso->lbmGetData( i   , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); 
 
140
                        *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); 
 
141
 
 
142
                        *this->mpIso->lbmGetData( i-1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); 
 
143
                        *this->mpIso->lbmGetData( i   , j   ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); 
 
144
                        *this->mpIso->lbmGetData( i+1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); 
 
145
 
 
146
                        *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); 
 
147
                        *this->mpIso->lbmGetData( i   , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); 
 
148
                        *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
111
149
        }
112
150
 
113
 
        
114
 
#if ELBEEM_PLUGIN!=1
115
 
        if(mUseTestdata) {
116
 
                int border = 1;
117
 
                for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++)
118
 
                        for(int j=0;j<mLevel[mMaxRefine].lSizey-1;j++) {
119
 
                                for(int l=0; l<=border; l++) {
120
 
                                        *D::mpIso->lbmGetData( l-1,                         j,ZKOFF) = *D::mpIso->lbmGetData( border+1,                         j,ZKOFF);
121
 
                                        *D::mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-l, j,ZKOFF) = *D::mpIso->lbmGetData( mLevel[mMaxRefine].lSizex-border-1, j,ZKOFF);
122
 
                                }
123
 
                        }
124
 
 
125
 
                for(int k=0;k<mLevel[mMaxRefine].lSizez-1;k++)
126
 
                        for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) {      
127
 
                                for(int l=0; l<=border; l++) {
128
 
                                        *D::mpIso->lbmGetData( i, l-1,                         ZKOFF) = *D::mpIso->lbmGetData( i, border+1,                         ZKOFF);
129
 
                                        *D::mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-l, ZKOFF) = *D::mpIso->lbmGetData( i, mLevel[mMaxRefine].lSizey-border-1, ZKOFF);
130
 
                                }
131
 
                        }
132
 
 
133
 
                if(D::cDimension == 3) {
134
 
                        // only for 3D
135
 
                        for(int j=-1;j<mLevel[mMaxRefine].lSizey+1;j++)
136
 
                                for(int i=-1;i<mLevel[mMaxRefine].lSizex+1;i++) {      
137
 
                                        for(int l=0; l<=border; l++) {
138
 
                                                *D::mpIso->lbmGetData( i,j,l-1                        ) = *D::mpIso->lbmGetData( i,j, border+1                         );
139
 
                                                *D::mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-l) = *D::mpIso->lbmGetData( i,j,mLevel[mMaxRefine].lSizez-1-border);
140
 
                                        }
141
 
                                }
142
 
                }
143
 
        } // testdata
144
 
#endif // ELBEEM_PLUGIN
145
 
        // */
146
151
 
147
152
        // update preview, remove 2d?
148
 
        if(D::mOutputSurfacePreview) {
149
 
                int pvsx = (int)(D::mPreviewFactor*D::mSizex);
150
 
                int pvsy = (int)(D::mPreviewFactor*D::mSizey);
151
 
                int pvsz = (int)(D::mPreviewFactor*D::mSizez);
152
 
                //float scale = (float)D::mSizex / previewSize;
153
 
                LbmFloat scalex = (LbmFloat)D::mSizex/(LbmFloat)pvsx;
154
 
                LbmFloat scaley = (LbmFloat)D::mSizey/(LbmFloat)pvsy;
155
 
                LbmFloat scalez = (LbmFloat)D::mSizez/(LbmFloat)pvsz;
156
 
                for(int k= 0; k< ((D::cDimension==3) ? (pvsz-1):1) ; ++k) 
 
153
        if((this->mOutputSurfacePreview)&&(LBMDIM==3)) {
 
154
                int pvsx = (int)(this->mPreviewFactor*this->mSizex);
 
155
                int pvsy = (int)(this->mPreviewFactor*this->mSizey);
 
156
                int pvsz = (int)(this->mPreviewFactor*this->mSizez);
 
157
                //float scale = (float)this->mSizex / previewSize;
 
158
                LbmFloat scalex = (LbmFloat)this->mSizex/(LbmFloat)pvsx;
 
159
                LbmFloat scaley = (LbmFloat)this->mSizey/(LbmFloat)pvsy;
 
160
                LbmFloat scalez = (LbmFloat)this->mSizez/(LbmFloat)pvsz;
 
161
                for(int k= 0; k< (pvsz-1); ++k) 
157
162
                for(int j=0;j< pvsy;j++) 
158
163
                for(int i=0;i< pvsx;i++) {
159
 
                                        *mpPreviewSurface->lbmGetData(i,j,k) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) );
 
164
                                        *mpPreviewSurface->lbmGetData(i,j,k) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) );
160
165
                                }
161
166
                // set borders again...
162
 
                for(int k= 0; k< ((D::cDimension == 3) ? (pvsz-1):1) ; ++k) {
 
167
                for(int k= 0; k< (pvsz-1); ++k) {
163
168
                        for(int j=0;j< pvsy;j++) {
164
 
                                *mpPreviewSurface->lbmGetData(0,j,k) = *D::mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) );
165
 
                                *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *D::mpIso->lbmGetData( D::mSizex-1, (int)(j*scaley), (int)(k*scalez) );
 
169
                                *mpPreviewSurface->lbmGetData(0,j,k) = *this->mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) );
 
170
                                *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *this->mpIso->lbmGetData( this->mSizex-1, (int)(j*scaley), (int)(k*scalez) );
166
171
                        }
167
172
                        for(int i=0;i< pvsx;i++) {
168
 
                                *mpPreviewSurface->lbmGetData(i,0,k) = *D::mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) );
169
 
                                *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *D::mpIso->lbmGetData( (int)(i*scalex), D::mSizey-1, (int)(k*scalez) );
 
173
                                *mpPreviewSurface->lbmGetData(i,0,k) = *this->mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) );
 
174
                                *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *this->mpIso->lbmGetData( (int)(i*scalex), this->mSizey-1, (int)(k*scalez) );
170
175
                        }
171
176
                }
172
 
                if(D::cDimension == 3) {
173
 
                        // only for 3D
 
177
                for(int j=0;j<pvsy;j++)
 
178
                        for(int i=0;i<pvsx;i++) {      
 
179
                                *mpPreviewSurface->lbmGetData(i,j,0) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0);
 
180
                                *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , this->mSizez-1);
 
181
                        }
 
182
 
 
183
                if(mFarFieldSize>=2.) {
 
184
                        // also remove preview border
 
185
                        for(int k= 0; k< (pvsz-1); ++k) {
 
186
                                for(int j=0;j< pvsy;j++) {
 
187
                                        *mpPreviewSurface->lbmGetData(0,j,k) = 
 
188
                                        *mpPreviewSurface->lbmGetData(1,j,k) =  
 
189
                                        *mpPreviewSurface->lbmGetData(2,j,k);
 
190
                                        *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = 
 
191
                                        *mpPreviewSurface->lbmGetData(pvsx-2,j,k) = 
 
192
                                        *mpPreviewSurface->lbmGetData(pvsx-3,j,k);
 
193
                                        //0.0;
 
194
                                }
 
195
                                for(int i=0;i< pvsx;i++) {
 
196
                                        *mpPreviewSurface->lbmGetData(i,0,k) = 
 
197
                                        *mpPreviewSurface->lbmGetData(i,1,k) = 
 
198
                                        *mpPreviewSurface->lbmGetData(i,2,k);
 
199
                                        *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = 
 
200
                                        *mpPreviewSurface->lbmGetData(i,pvsy-2,k) = 
 
201
                                        *mpPreviewSurface->lbmGetData(i,pvsy-3,k);
 
202
                                        //0.0;
 
203
                                }
 
204
                        }
174
205
                        for(int j=0;j<pvsy;j++)
175
206
                                for(int i=0;i<pvsx;i++) {      
176
 
                                        *mpPreviewSurface->lbmGetData(i,j,0) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0);
177
 
                                        *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *D::mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , D::mSizez-1);
178
 
                                }
179
 
                } // borders done...
 
207
                                        *mpPreviewSurface->lbmGetData(i,j,0) = 
 
208
                                        *mpPreviewSurface->lbmGetData(i,j,1) = 
 
209
                                        *mpPreviewSurface->lbmGetData(i,j,2);
 
210
                                        *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = 
 
211
                                        *mpPreviewSurface->lbmGetData(i,j,pvsz-2) = 
 
212
                                        *mpPreviewSurface->lbmGetData(i,j,pvsz-3);
 
213
                                        //0.0;
 
214
                        }
 
215
                }
180
216
        }
181
217
 
182
 
#if ELBEEM_PLUGIN!=1
183
 
        if(D::mInitDone) { 
184
 
                if(mpTest->mDebugvalue3<=0.0) handleTestdata();
185
 
        }
186
 
#endif // ELBEEM_PLUGIN!=1
187
218
        // correction
188
219
        return;
189
220
}
190
221
 
191
222
/*! calculate speeds of fluid objects (or inflow) */
192
 
template<class D>
193
 
void LbmFsgrSolver<D>::recalculateObjectSpeeds() {
194
 
        int numobjs = (int)(D::mpGiObjects->size());
 
223
void LbmFsgrSolver::recalculateObjectSpeeds() {
 
224
        const bool debugRecalc = false;
 
225
        int numobjs = (int)(this->mpGiObjects->size());
195
226
        // note - (numobjs + 1) is entry for domain settings
 
227
 
 
228
        if(debugRecalc) errMsg("recalculateObjectSpeeds","start, #obj:"<<numobjs);
196
229
        if(numobjs>255-1) {
197
230
                errFatal("LbmFsgrSolver::recalculateObjectSpeeds","More than 256 objects currently not supported...",SIMWORLD_INITERROR);
198
231
                return;
199
232
        }
200
233
        mObjectSpeeds.resize(numobjs+1);
201
234
        for(int i=0; i<(int)(numobjs+0); i++) {
202
 
                mObjectSpeeds[i] = vec2L(D::mpParam->calculateLattVelocityFromRw( vec2P( (*D::mpGiObjects)[i]->getInitialVelocity() )));
203
 
                //errMsg("recalculateObjectSpeeds","id"<<i<<" set to "<< mObjectSpeeds[i]<<", unscaled:"<< (*D::mpGiObjects)[i]->getInitialVelocity() );
 
235
                mObjectSpeeds[i] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) )));
 
236
                if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" set to "<< mObjectSpeeds[i]<<", unscaled:"<< (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) );
204
237
        }
205
238
 
206
239
        // also reinit part slip values here
207
240
        mObjectPartslips.resize(numobjs+1);
208
 
        for(int i=0; i<(int)(numobjs+0); i++) {
209
 
                mObjectPartslips[i] = (LbmFloat)(*D::mpGiObjects)[i]->getGeoPartSlipValue();
 
241
        for(int i=0; i<=(int)(numobjs+0); i++) {
 
242
                if(i<numobjs) {
 
243
                        mObjectPartslips[i] = (LbmFloat)(*this->mpGiObjects)[i]->getGeoPartSlipValue();
 
244
                } else {
 
245
                        // domain setting
 
246
                        mObjectPartslips[i] = this->mDomainPartSlipValue;
 
247
                }
 
248
                LbmFloat set = mObjectPartslips[i];
 
249
 
 
250
                // as in setInfluenceVelocity
 
251
                const LbmFloat dt = mLevel[mMaxRefine].timestep;
 
252
                const LbmFloat dtInter = 0.01;
 
253
                //LbmFloat facFv = 1.-set;
 
254
                // mLevel[mMaxRefine].timestep
 
255
                LbmFloat facNv = (LbmFloat)( 1.-pow( (double)(set), (double)(dt/dtInter)) );
 
256
                errMsg("mObjectPartslips","id:"<<i<<"/"<<numobjs<<"  ts:"<<dt<< " its:"<<(dt/dtInter) <<" set"<<set<<" nv"<<facNv<<" test:"<<
 
257
                         pow( (double)(1.-facNv),(double)(dtInter/dt))  );
 
258
                mObjectPartslips[i] = facNv;
 
259
 
 
260
                if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" parts "<< mObjectPartslips[i] );
210
261
        }
211
 
        //errMsg("GEOIN"," dm set "<<mDomainPartSlipValue);
212
 
        mObjectPartslips[numobjs] = mDomainPartSlipValue;
 
262
 
 
263
        if(debugRecalc) errMsg("recalculateObjectSpeeds","done, domain:"<<mObjectPartslips[numobjs]<<" n"<<numobjs);
213
264
}
214
265
 
215
266
 
217
268
/*****************************************************************************/
218
269
/*! debug object display */
219
270
/*****************************************************************************/
220
 
template<class D>
221
 
vector<ntlGeometryObject*> LbmFsgrSolver<D>::getDebugObjects() { 
 
271
vector<ntlGeometryObject*> LbmFsgrSolver::getDebugObjects() { 
222
272
        vector<ntlGeometryObject*> debo; 
223
 
        if(D::mOutputSurfacePreview) {
 
273
        if(this->mOutputSurfacePreview) {
224
274
                debo.push_back( mpPreviewSurface );
225
275
        }
226
 
#if ELBEEM_PLUGIN!=1
 
276
#if LBM_INCLUDE_TESTSOLVERS==1
227
277
        if(mUseTestdata) {
228
278
                vector<ntlGeometryObject*> tdebo; 
229
279
                tdebo = mpTest->getDebugObjects();
238
288
 *****************************************************************************/
239
289
 
240
290
/*! init particle positions */
241
 
template<class D>
242
 
int LbmFsgrSolver<D>::initParticles(ParticleTracer *partt) { 
243
 
#if ELBEEM_PLUGIN==1
244
 
        partt = NULL; // remove warning
245
 
#else // ELBEEM_PLUGIN
 
291
int LbmFsgrSolver::initParticles() { 
246
292
  int workSet = mLevel[mMaxRefine].setCurr;
247
293
  int tries = 0;
248
294
  int num = 0;
249
 
        mpParticles=partt;
250
 
 
251
 
  //partt->setSimEnd  ( ntlVec3Gfx(D::mSizex-1, D::mSizey-1, getForZMax1()) );
252
 
  partt->setSimEnd  ( ntlVec3Gfx(D::mSizex,   D::mSizey,   getForZMaxBnd(mMaxRefine)) );
 
295
        ParticleTracer *partt = mpParticles;
 
296
 
 
297
  partt->setStart( this->mvGeoStart + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
 
298
  partt->setEnd  ( this->mvGeoEnd   + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
 
299
 
253
300
  partt->setSimStart( ntlVec3Gfx(0.0) );
 
301
  partt->setSimEnd  ( ntlVec3Gfx(this->mSizex,   this->mSizey,   getForZMaxBnd(mMaxRefine)) );
254
302
  
255
 
  while( (num<partt->getNumParticles()) && (tries<100*partt->getNumParticles()) ) {
256
 
    double x,y,z;
257
 
    x = 0.0+(( (float)(D::mSizex-1) )     * (rand()/(RAND_MAX+1.0)) );
258
 
    y = 0.0+(( (float)(D::mSizey-1) )     * (rand()/(RAND_MAX+1.0)) );
259
 
    z = 0.0+(( (float) getForZMax1(mMaxRefine) )* (rand()/(RAND_MAX+1.0)) );
260
 
    int i = (int)(x-0.5);
261
 
    int j = (int)(y-0.5);
262
 
    int k = (int)(z-0.5);
263
 
    if(D::cDimension==2) {
264
 
      k = 0;
265
 
      z = 0.5; // place in the middle of domain
 
303
  while( (num<partt->getNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) {
 
304
    LbmFloat x,y,z,t;
 
305
    x = 1.0+(( (LbmFloat)(this->mSizex-3.) )          * (rand()/(RAND_MAX+1.0)) );
 
306
    y = 1.0+(( (LbmFloat)(this->mSizey-3.) )          * (rand()/(RAND_MAX+1.0)) );
 
307
    z = 1.0+(( (LbmFloat) getForZMax1(mMaxRefine)-2. )* (rand()/(RAND_MAX+1.0)) );
 
308
    int i = (int)(x+0.5);
 
309
    int j = (int)(y+0.5);
 
310
    int k = (int)(z+0.5);
 
311
    if(LBMDIM==2) {
 
312
      k = 0; z = 0.5; // place in the middle of domain
266
313
    }
267
314
 
268
 
    if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ||
269
 
        TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) { // only fluid cells?
 
315
    //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid) ) 
 
316
    //&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid ) 
 
317
    //if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) { 
 
318
    if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid) ) { 
 
319
                        bool cellOk = true;
 
320
                        //? FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; }
 
321
                        if(!cellOk) continue;
270
322
      // in fluid...
271
323
      partt->addParticle(x,y,z);
 
324
                        partt->getLast()->setStatus(PART_IN);
 
325
                        partt->getLast()->setType(PART_TRACER);
 
326
 
 
327
                        partt->getLast()->setSize(1.);
 
328
                        // randomize size
 
329
                        partt->getLast()->setSize(0.5 + (rand()/(RAND_MAX+1.0)));
 
330
 
 
331
                        if( ( partt->getInitStart()>0.)
 
332
                                        && ( partt->getInitEnd()>0.)
 
333
                                        && ( partt->getInitEnd()>partt->getInitStart() )) {
 
334
                t = partt->getInitStart()+ (partt->getInitEnd()-partt->getInitStart())*(rand()/(RAND_MAX+1.0));
 
335
                                partt->getLast()->setLifeTime( -t );
 
336
                        }
272
337
      num++;
273
338
    }
274
339
    tries++;
275
 
  }
276
 
  debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles ", 10);
 
340
  } // */
 
341
 
 
342
        /*FSGR_FORIJK1(mMaxRefine) {
 
343
                if( (RFLAG(mMaxRefine,i,j,k, workSet) & (CFNoBndFluid)) ) {
 
344
        LbmFloat rndn = (rand()/(RAND_MAX+1.0));
 
345
                        if(rndn>0.0) {
 
346
                                ntlVec3Gfx pos( (LbmFloat)(i)-0.5, (LbmFloat)(j)-0.5, (LbmFloat)(k)-0.5 );
 
347
                                if(LBMDIM==2) { pos[2]=0.5; }
 
348
                                partt->addParticle( pos[0],pos[1],pos[2] );
 
349
                                partt->getLast()->setStatus(PART_IN);
 
350
                                partt->getLast()->setType(PART_TRACER);
 
351
                                partt->getLast()->setSize(1.0);
 
352
                        }
 
353
                }
 
354
        } // */
 
355
 
 
356
 
 
357
        // DEBUG TEST
 
358
#if LBM_INCLUDE_TESTSOLVERS==1
 
359
        if(mUseTestdata) { 
 
360
                const bool partDebug=false;
 
361
                if(mpTest->mDebugvalue2!=0.0){ errMsg("LbmTestdata"," part init "<<mpTest->mDebugvalue2); }
 
362
                if(mpTest->mDebugvalue2==-12.0){ 
 
363
                        const int lev = mMaxRefine;
 
364
                        for(int i=5;i<15;i++) {
 
365
                                LbmFloat x,y,z;
 
366
                                y = 0.5+(LbmFloat)(i);
 
367
                                x = mLevel[lev].lSizex/20.0*10.0;
 
368
                                z = mLevel[lev].lSizez/20.0*2.0;
 
369
                                partt->addParticle(x,y,z);
 
370
                                partt->getLast()->setStatus(PART_IN);
 
371
                                partt->getLast()->setType(PART_BUBBLE);
 
372
                                partt->getLast()->setSize(  (-4.0+(LbmFloat)i)/1.0  );
 
373
                                if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
 
374
                        }
 
375
                }
 
376
                if(mpTest->mDebugvalue2==-11.0){ 
 
377
                        const int lev = mMaxRefine;
 
378
                        for(int i=5;i<15;i++) {
 
379
                                LbmFloat x,y,z;
 
380
                                y = 10.5+(LbmFloat)(i);
 
381
                                x = mLevel[lev].lSizex/20.0*10.0;
 
382
                                z = mLevel[lev].lSizez/20.0*40.0;
 
383
                                partt->addParticle(x,y,z);
 
384
                                partt->getLast()->setStatus(PART_IN);
 
385
                                partt->getLast()->setType(PART_DROP);
 
386
                                partt->getLast()->setSize(  (-4.0+(LbmFloat)i)/1.0  );
 
387
                                if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
 
388
                        }
 
389
                }
 
390
                // place floats on rectangular region FLOAT_JITTER_BND
 
391
                if(mpTest->mDebugvalue2==-10.0){ 
 
392
                        const int lev = mMaxRefine;
 
393
                        const int sx = mLevel[lev].lSizex;
 
394
                        const int sy = mLevel[lev].lSizey;
 
395
                        //for(int j=-(int)(sy*0.25);j<-(int)(sy*0.25)+2;++j) { for(int i=-(int)(sx*0.25);i<-(int)(sy*0.25)+2;++i) {
 
396
                        //for(int j=-(int)(sy*1.25);j<(int)(2.25*sy);++j) { for(int i=-(int)(sx*1.25);i<(int)(2.25*sx);++i) {
 
397
                        for(int j=-(int)(sy*0.3);j<(int)(1.3*sy);++j) { for(int i=-(int)(sx*0.3);i<(int)(1.3*sx);++i) {
 
398
                        //for(int j=-(int)(sy*0.2);j<(int)(0.2*sy);++j) { for(int i= (int)(sx*0.5);i<= (int)(0.51*sx);++i) {
 
399
                                        LbmFloat x,y,z;
 
400
                                        x = 0.0+(LbmFloat)(i);
 
401
                                        y = 0.0+(LbmFloat)(j);
 
402
                                        //z = mLevel[lev].lSizez/10.0*2.5 - 1.0;
 
403
                                        z = mLevel[lev].lSizez/20.0*9.5 - 1.0;
 
404
                                        //z = mLevel[lev].lSizez/20.0*4.5 - 1.0;
 
405
                                        partt->addParticle(x,y,z);
 
406
                                        //if( (i>0)&&(i<sx) && (j>0)&&(j<sy) ) { partt->getLast()->setStatus(PART_IN); } else { partt->getLast()->setStatus(PART_OUT); }
 
407
                                        partt->getLast()->setStatus(PART_IN);
 
408
                                        partt->getLast()->setType(PART_FLOAT);
 
409
                                        partt->getLast()->setSize( 15.0 );
 
410
                                        if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
 
411
                         }
 
412
                }       }
 
413
        } 
 
414
        // DEBUG TEST
 
415
#endif // LBM_INCLUDE_TESTSOLVERS
 
416
 
 
417
        
 
418
  debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb<<", tries:"<<tries, 10);
277
419
  if(num != partt->getNumParticles()) return 1;
278
 
#endif // ELBEEM_PLUGIN
279
420
 
280
421
        return 0;
281
422
}
282
423
 
283
 
template<class D>
284
 
void LbmFsgrSolver<D>::advanceParticles(ParticleTracer *partt) { 
285
 
#if ELBEEM_PLUGIN==1
286
 
        partt = NULL; // remove warning
287
 
#else // ELBEEM_PLUGIN
 
424
#define P_CHANGETYPE(p, newtype) \
 
425
                p->setLifeTime(0.); \
 
426
    /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \
 
427
                p->setType(newtype); 
 
428
 
 
429
void LbmFsgrSolver::advanceParticles() { 
288
430
  int workSet = mLevel[mMaxRefine].setCurr;
289
431
        LbmFloat vx=0.0,vy=0.0,vz=0.0;
290
432
        LbmFloat rho, df[27]; //feq[27];
291
 
        if(mpParticles!=partt) { errMsg("LbmFsgrSolver<D>::advanceParticles","Invalid ParticleTracer..."); }
 
433
#define DEL_PART { \
 
434
        /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<"  ");  */ \
 
435
        p->setActive( false ); \
 
436
        continue; }
 
437
 
 
438
        myTime_t parttstart = getTime(); 
 
439
        const LbmFloat cellsize = this->mpParam->getCellSize();
 
440
        const LbmFloat timestep = this->mpParam->getTimestep();
 
441
        //const LbmFloat viscAir = 1.79 * 1e-5; // RW2L kin. viscosity, mu
 
442
        const LbmFloat viscWater = 1.0 * 1e-6; // RW2L kin. viscosity, mu
 
443
        const LbmFloat rhoAir = 1.2;  // [kg m^-3] RW2L
 
444
        const LbmFloat rhoWater = 1000.0; // RW2L
 
445
        const LbmFloat minDropSize = 0.0005; // [m], = 2mm  RW2L
 
446
        const LbmVec   velAir(0.); // [m / s]
 
447
 
 
448
        const LbmFloat r1 = 0.005;  // r max
 
449
        const LbmFloat r2 = 0.0005; // r min
 
450
        const LbmFloat v1 = 9.0; // v max
 
451
        const LbmFloat v2 = 2.0; // v min
 
452
        const LbmVec rwgrav = vec2L( this->mpParam->getGravity(mSimulationTime) );
 
453
        const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){ 
292
454
 
293
455
        // TODO use timestep size
294
 
  for(vector<ParticleObject>::iterator pit= partt->getParticlesBegin();
295
 
      pit!= partt->getParticlesEnd(); pit++) {
296
 
    //errorOut(" pit "<< (*pit).getPos() );
 
456
        //bool isIn,isOut,isInZ;
 
457
        //const int cutval = 1+mCutoff/2; // TODO FIXME add half border!
 
458
        //const int cutval = mCutoff/2; // TODO FIXME add half border!
 
459
        //const int cutval = 0; // TODO FIXME add half border!
 
460
        const int cutval = mCutoff; // use full border!?
 
461
        int actCnt=0;
 
462
        if(this->mStepCnt%50==49) { mpParticles->cleanup(); }
 
463
  for(vector<ParticleObject>::iterator pit= mpParticles->getParticlesBegin();
 
464
      pit!= mpParticles->getParticlesEnd(); pit++) {
 
465
    //errMsg("PIT"," pit"<<(*pit).getId()<<" pos:"<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel()  );
 
466
    //errMsg("PIT"," pit pos:"<< (*pit)->getPos()<<" vel:"<< (*pit)->getVel()<<" status:"<<convertFlags2String((*pit)->getFlags()) <<" " <<mpParticles->getStart()<<" "<<mpParticles->getEnd() );
 
467
                //int flag = (*pit).getFlags();
297
468
    if( (*pit).getActive()==false ) continue;
 
469
                // skip until reached
 
470
                ParticleObject *p = &(*pit);
 
471
                if(p->getLifeTime()<0.){ 
 
472
                        if(p->getLifeTime()<-this->mSimulationTime) continue; 
 
473
                        else p->setLifeTime(-mLevel[mMaxRefine].timestep); // zero for following update
 
474
                }
298
475
    int i,j,k;
299
 
                ParticleObject *p = &(*pit);
 
476
                p->setLifeTime(p->getLifeTime()+mLevel[mMaxRefine].timestep);
300
477
 
301
478
                // nearest neighbor, particle positions don't include empty bounds
302
479
                ntlVec3Gfx pos = p->getPos();
303
480
                i= (int)(pos[0]+0.5);
304
481
                j= (int)(pos[1]+0.5);
305
482
                k= (int)(pos[2]+0.5);
306
 
                if(D::cDimension==2) {
307
 
                        k = 0;
308
 
                }
309
 
 
310
 
                if( (i<0)||(i>D::mSizex-1)||
311
 
                                (j<0)||(j>D::mSizey-1)||
312
 
                                (k<0)||(k>D::mSizez-1) ) {
313
 
                        p->setActive( false );
314
 
                        continue;
315
 
                }
316
 
 
317
 
                if(p->getStatus()==0) {
 
483
                if(LBMDIM==2) { k = 0; }
 
484
 
 
485
                // only testdata handling, all for sws
 
486
#if LBM_INCLUDE_TESTSOLVERS==1
 
487
                if(useff && (mpTest->mDebugvalue1>0.)) {
 
488
                        p->setStatus(PART_OUT);
 
489
                        mpTest->handleParticle(p, i,j,k); continue;
 
490
                } 
 
491
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
492
 
 
493
                // FIXME , add k tests again, remove per type ones...
 
494
                if(p->getStatus()&PART_IN) { // IN
 
495
                        if( (i<cutval)||(i>this->mSizex-1-cutval)||
 
496
                                        (j<cutval)||(j>this->mSizey-1-cutval)
 
497
                                        //||(k<cutval)||(k>this->mSizez-1-cutval) 
 
498
                                        ) {
 
499
                                if(!useff) { DEL_PART;
 
500
                                } else { 
 
501
                                        p->setStatus(PART_OUT); 
 
502
                                        /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART;
 
503
                                }
 
504
                        } 
 
505
                } else { // OUT rough check
 
506
                        // check in again?
 
507
                        if( (i>=cutval)&&(i<=this->mSizex-1-cutval)&&
 
508
                                        (j>=cutval)&&(j<=this->mSizey-1-cutval)
 
509
                                        //&&(k>=cutval)&&(k<=this->mSizez-1-cutval) 
 
510
                                        ) {
 
511
                                p->setStatus(PART_IN);
 
512
                                /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART;
 
513
                        }
 
514
                }
 
515
                //p->setStatus(PART_OUT);// DEBUG always out!
 
516
 
 
517
                if( (p->getType()==PART_BUBBLE) ||
 
518
                    (p->getType()==PART_TRACER) ) {
 
519
 
318
520
                        // no interpol
319
521
                        rho = vx = vy = vz = 0.0;
320
 
                        FORDF0{
321
 
                                LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
322
 
                                df[l] = cdf;
323
 
                                rho += cdf; 
324
 
                                vx  += (D::dfDvecX[l]*cdf); 
325
 
                                vy  += (D::dfDvecY[l]*cdf);  
326
 
                                vz  += (D::dfDvecZ[l]*cdf);  
327
 
                        }
328
 
 
329
 
                        // remove gravity influence
330
 
                        //FORDF0{ feq[l] = D::getCollideEq(l, rho,vx,vy,vz); }
331
 
                        //const LbmFloat Qo = D::getLesNoneqTensorCoeff(df,feq);
332
 
                        //const LbmFloat lesomega = D::getLesOmega(mLevel[mMaxRefine].omega,mLevel[mMaxRefine].lcsmago,Qo);
333
 
                        const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
334
 
                        vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
335
 
                        vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
336
 
                        vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
337
 
 
338
 
                        if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ||
339
 
                                        TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
340
 
                                // still ok
341
 
                        } else {
342
 
                                // out of bounds, deactivate...
343
 
                                // FIXME make fsgr treatment
344
 
                                p->setActive( false );
345
 
                                continue;
346
 
                                D::mNumParticlesLost++;
347
 
                        }
348
 
 
349
 
                        p->advance( vx,vy,vz );
350
 
                  // fluid particle
351
 
                } else {
352
 
                        p->setVel( p->getVel() * 0.999 ); // dampen...
353
 
                        p->addToVel( vec2G(mLevel[mMaxRefine].gravity) );
354
 
                        p->advanceVel();
355
 
                        //errMsg("NNNPART"," at "<<p->getPos()<<" u="<<p->getVel() );
356
 
                        if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ||
357
 
                                        TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
358
 
                                // still ok
359
 
                        } else {
360
 
                                // out of bounds, deactivate...
361
 
                                // FIXME make fsgr treatment
362
 
                                p->setActive( false );
363
 
                                continue;
364
 
                                D::mNumParticlesLost++;
365
 
                        }
 
522
                        if(p->getStatus()&PART_IN) { // IN
 
523
                                if(k>=cutval) {
 
524
                                        if(k>this->mSizez-1-cutval) DEL_PART; 
 
525
                                        FORDF0{
 
526
                                                LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
 
527
                                                df[l] = cdf;
 
528
                                                rho += cdf; 
 
529
                                                vx  += (this->dfDvecX[l]*cdf); 
 
530
                                                vy  += (this->dfDvecY[l]*cdf);  
 
531
                                                vz  += (this->dfDvecZ[l]*cdf);  
 
532
                                        }
 
533
 
 
534
                                        // remove gravity influence
 
535
                                        const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
 
536
                                        vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
 
537
                                        vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
 
538
                                        vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
 
539
 
 
540
                                        if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid) ) {
 
541
                                                // still ok
 
542
                                        } else { // OUT
 
543
                                                // out of bounds, deactivate...
 
544
                                                // FIXME make fsgr treatment
 
545
                                                if(p->getType()==PART_BUBBLE) { P_CHANGETYPE(p, PART_FLOAT ); continue; }
 
546
                                        }
 
547
                                } else {
 
548
                                        // below 3d region, just rise
 
549
                                }
 
550
                        } else { // OUT
 
551
#if LBM_INCLUDE_TESTSOLVERS==1
 
552
                                if(useff) { mpTest->handleParticle(p, i,j,k); }
 
553
                                else DEL_PART;
 
554
#else // LBM_INCLUDE_TESTSOLVERS==1
 
555
                                DEL_PART;
 
556
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
557
                                // TODO use x,y vel...?
 
558
                        }
 
559
 
 
560
                        ntlVec3Gfx v = p->getVel(); // dampen...
 
561
                        if( (useff)&& (p->getType()==PART_BUBBLE) ) {
 
562
                                // test rise
 
563
                                //O vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
 
564
 
 
565
                                LbmFloat radius = p->getSize() * minDropSize;
 
566
                                LbmVec   velPart = vec2L(p->getVel()) *cellsize/timestep; // L2RW, lattice velocity
 
567
                                LbmVec   velWater = LbmVec(vx,vy,vz) *cellsize/timestep;// L2RW, fluid velocity
 
568
                                LbmVec   velRel = velWater - velPart;
 
569
                                LbmFloat velRelNorm = norm(velRel);
 
570
                                // TODO calculate values in lattice units, compute CD?!??!
 
571
                                LbmFloat pvolume = rhoAir * 4.0/3.0 * M_PI* radius*radius*radius; // volume: 4/3 pi r^3
 
572
                                //const LbmFloat cd = 
 
573
 
 
574
                                LbmVec fb = -rwgrav* pvolume *rhoWater;
 
575
                                LbmVec fd = velRel*6.0*M_PI*radius* (1e-3); //viscWater;
 
576
                                LbmVec change = (fb+fd) *10.0*timestep  *(timestep/cellsize);
 
577
                                //LbmVec change = (fb+fd) *timestep / (pvolume*rhoAir)  *(timestep/cellsize);
 
578
                                //actCnt++; // should be after active test
 
579
                                if(actCnt<0) {
 
580
                                        errMsg("\nPIT","BTEST1   vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
 
581
                                        errMsg("PIT","BTEST2        cellsize="<<cellsize<<" timestep="<<timestep<<" viscW="<<viscWater<<" ss/mb="<<(timestep/(pvolume*rhoAir)));
 
582
                                        errMsg("PIT","BTEST2        grav="<<rwgrav<<"  " );
 
583
                                        errMsg("PIT","BTEST2        change="<<(change)<<" fb="<<(fb)<<" fd="<<(fd)<<" ");
 
584
                                        errMsg("PIT","BTEST2        change="<<norm(change)<<" fb="<<norm(fb)<<" fd="<<norm(fd)<<" ");
 
585
#if LOOPTEST==1
 
586
                                        errMsg("PIT","BTEST2        n="<<n<<" "); // LOOPTEST! DEBUG
 
587
#endif // LOOPTEST==1
 
588
                                        errMsg("PIT","\n");
 
589
                                }
 
590
                                        
 
591
                                //v += change;
 
592
                                //v += ntlVec3Gfx(vx,vy,vz);
 
593
                                LbmVec fd2 = (LbmVec(vx,vy,vz)-vec2L(p->getVel())) * 6.0*M_PI*radius* (1e-3); //viscWater;
 
594
                                LbmFloat w = 0.99;
 
595
                                vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[mMaxRefine].gravity[2]);
 
596
                                v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2);
 
597
                        } else if(p->getType()==PART_TRACER) {
 
598
                                v = ntlVec3Gfx(vx,vy,vz);
 
599
                                CellFlagType fflag = RFLAG(mMaxRefine, i,j,k, workSet);
 
600
 
 
601
                                if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true);
 
602
                                } else { p->setInFluid(false); }
 
603
 
 
604
                                //if(( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) {
 
605
                                if( (( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) ||
 
606
                                                (( fflag&CFInter ) && (!(fflag&CFNoNbFluid)) ) ) {
 
607
                                        // only real fluid
 
608
#define TRACE_JITTER 0.025
 
609
#define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
 
610
#define __TRACE_RAND 0.
 
611
#if LBMDIM==3
 
612
                                        p->advance( TRACE_RAND,TRACE_RAND,TRACE_RAND);
 
613
#else
 
614
                                        p->advance( TRACE_RAND,TRACE_RAND, 0.);
 
615
#endif
 
616
 
 
617
                                //} // test!?
 
618
                                //if( fflag&(CFNoBndFluid) ) {
 
619
                                        // ok
 
620
                                //} else if( fflag&(CFEmpty|CFNoNbFluid) ) {
 
621
                                //} else if( fflag&(CFEmpty) ) {
 
622
                                        //v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
 
623
                                } else {
 
624
                                        // move inwards along normal, make sure normal is valid first
 
625
                                        const int lev = mMaxRefine;
 
626
                                        LbmFloat nx=0.,ny=0.,nz=0., nv1,nv2;
 
627
                                        bool nonorm = false;
 
628
                                        if(i<=0)              { nx = -1.; nonorm = true; }
 
629
                                        if(i>=this->mSizex-1) { nx =  1.; nonorm = true; }
 
630
                                        if(j<=0)              { ny = -1.; nonorm = true; }
 
631
                                        if(j>=this->mSizey-1) { ny =  1.; nonorm = true; }
 
632
#if LBMDIM==3
 
633
                                        if(k<=0)              { nz = -1.; nonorm = true; }
 
634
                                        if(k>=this->mSizez-1) { nz =  1.; nonorm = true; }
 
635
#endif // LBMDIM==3
 
636
                        //mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
 
637
                                        if(!nonorm) {
 
638
                                                if(RFLAG_NB(lev,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dE,dFfrac); } else nv1 = 0.0;
 
639
                                                if(RFLAG_NB(lev,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dW,dFfrac); } else nv2 = 0.0;
 
640
                                                nx = 0.5* (nv2-nv1);
 
641
                                                if(RFLAG_NB(lev,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dN,dFfrac); } else nv1 = 0.0;
 
642
                                                if(RFLAG_NB(lev,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dS,dFfrac); } else nv2 = 0.0;
 
643
                                                ny = 0.5* (nv2-nv1);
 
644
#if LBMDIM==3
 
645
                                                if(RFLAG_NB(lev,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dT,dFfrac); } else nv1 = 0.0;
 
646
                                                if(RFLAG_NB(lev,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dB,dFfrac); } else nv2 = 0.0;
 
647
                                                nz = 0.5* (nv2-nv1);
 
648
#else //LBMDIM==3
 
649
                                                nz = 0.0;
 
650
#endif //LBMDIM==3
 
651
                                        } else {
 
652
                                                v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
 
653
                                        }
 
654
 
 
655
                                        p->advanceVec( (ntlVec3Gfx(nx,ny,nz)) * -0.1 ); // + vec2G(mLevel[mMaxRefine].gravity);
 
656
                                        //if(normNoSqrt(v)<LBM_EPSILON) v = mLevel[mMaxRefine].gravity;
 
657
                                }
 
658
                        }
 
659
 
 
660
                        p->setVel( v );
 
661
                        p->advanceVel();
 
662
                        //errMsg("PPPP"," pos"<<p->getPos()<<" "<<p->getVel() );
 
663
                } 
 
664
 
 
665
                // drop handling
 
666
                else if(p->getType()==PART_DROP) {
 
667
                        ntlVec3Gfx v = p->getVel(); // dampen...
 
668
 
 
669
                        if(1) {
 
670
                                LbmFloat radius = p->getSize() * minDropSize;
 
671
                                LbmVec   velPart = vec2L(p->getVel()) *cellsize /timestep; // * cellsize / timestep; // L2RW, lattice velocity
 
672
                                LbmVec   velRel = velAir - velPart;
 
673
                                //LbmVec   velRelLat = velRel /cellsize*timestep; // L2RW
 
674
                                LbmFloat velRelNorm = norm(velRel);
 
675
                                // TODO calculate values in lattice units, compute CD?!??!
 
676
                                LbmFloat mb = rhoWater * 4.0/3.0 * M_PI* radius*radius*radius; // mass: 4/3 pi r^3 rho
 
677
                                const LbmFloat rw = (r1-radius)/(r1-r2);
 
678
                                const LbmFloat rmax = (0.5 + 0.5*rw);
 
679
                                const LbmFloat vmax = (v2 + (v1-v2)* (1.0-rw) );
 
680
                                const LbmFloat cd = (rmax) * (velRelNorm)/(vmax);
 
681
 
 
682
                                LbmVec fg = rwgrav * mb;//  * (1.0-rhoAir/rhoWater);
 
683
                                LbmVec fd = velRel* velRelNorm* cd*M_PI *rhoAir *0.5 *radius*radius;
 
684
                                LbmVec change = (fg+   fd ) *timestep / mb  *(timestep/cellsize);
 
685
 
 
686
                                //actCnt++; // should be after active test
 
687
                                if(actCnt<0) {
 
688
                                        errMsg("\nPIT","NTEST1   mb="<<mb<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel<<" pgetVel="<<p->getVel() );
 
689
                                        //errMsg("PIT","NTEST2        cellsize="<<cellsize<<" timestep="<<timestep<<" viscAir="<<viscAir<<" ss/mb="<<(timestep/mb));
 
690
                                        //errMsg("PIT","NTEST2        grav="<<rwgrav<<" mb="<<mb<<" "<<" cd="<<cd );
 
691
                                        //errMsg("PIT","NTEST2        change="<<norm(change)<<" fg="<<norm(fg)<<" fd="<<norm(fd)<<" ");
 
692
                                }
 
693
 
 
694
                                v += vec2G(change);
 
695
                                p->setVel(v); 
 
696
                                // NEW
 
697
                        } else {
 
698
                                p->setVel( v * 0.999 ); // dampen...
 
699
                                p->setVel( v ); // DEBUG!
 
700
                                p->addToVel( vec2G(mLevel[mMaxRefine].gravity) );\
 
701
                        } // OLD
 
702
                        p->advanceVel();
 
703
 
 
704
                        if(p->getStatus()&PART_IN) { // IN
 
705
                                if(k<cutval) { DEL_PART; continue; }
 
706
                                if(k<=this->mSizez-1-cutval){ 
 
707
                                        //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter)) {
 
708
                                        if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) {
 
709
                                                // still ok
 
710
                                        // shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
 
711
                                        } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid) ){
 
712
                                                // FIXME make fsgr treatment
 
713
                                                P_CHANGETYPE(p, PART_FLOAT ); continue; 
 
714
                                                // jitter in cell to prevent stacking when hitting a steep surface
 
715
                                                ntlVec3Gfx cpos = p->getPos(); cpos[0] += (rand()/(RAND_MAX+1.0))-0.5;
 
716
                                                cpos[1] += (rand()/(RAND_MAX+1.0))-0.5; p->setPos(cpos);
 
717
                                                //} else DEL_PART;
 
718
                                        } else {
 
719
                                                DEL_PART;
 
720
                                                this->mNumParticlesLost++;
 
721
                                        }
 
722
                                }
 
723
                        } else { // OUT
 
724
#if LBM_INCLUDE_TESTSOLVERS==1
 
725
                                if(useff) { mpTest->handleParticle(p, i,j,k); }
 
726
                                else{ DEL_PART; }
 
727
#else // LBM_INCLUDE_TESTSOLVERS==1
 
728
                                { DEL_PART; }
 
729
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
730
                        }
 
731
 
366
732
                } // air particle
 
733
 
 
734
                // inter particle
 
735
                else if(p->getType()==PART_INTER) {
 
736
                        if(p->getStatus()&PART_IN) { // IN
 
737
                                if((k<cutval)||(k>this->mSizez-1-cutval)) {
 
738
                                        // undecided particle above or below... remove?
 
739
                                        DEL_PART; 
 
740
                                }
 
741
 
 
742
                                if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
 
743
                                        // still ok
 
744
                                } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
 
745
                        //errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
746
 
 
747
                                        //P_CHANGETYPE(p, PART_BUBBLE ); continue;
 
748
                                        // currently bubbles off! DEBUG!
 
749
                                        //DEL_PART; // DEBUG bubbles off for siggraph
 
750
                                        P_CHANGETYPE(p, PART_FLOAT ); continue;
 
751
                                } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
 
752
                        //errMsg("PIT","NEWDROP pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
753
                                        P_CHANGETYPE(p, PART_DROP ); continue;
 
754
                                } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFBnd ) ) {
 
755
                                        // 060423 new test
 
756
                                        P_CHANGETYPE(p, PART_FLOAT ); continue;
 
757
                                }
 
758
                        } else { // OUT
 
759
                                // undecided particle outside... remove?
 
760
                                DEL_PART; 
 
761
                                //P_CHANGETYPE(p, PART_FLOAT ); continue;
 
762
                        }
 
763
                }
 
764
 
 
765
                // float particle
 
766
                else if(p->getType()==PART_FLOAT) {
 
767
                        //  test - delte on boundary!?
 
768
                        //if( (i<=cutval)||(i>=this->mSizex-1-cutval)|| (j<=cutval)||(j>=this->mSizey-1-cutval)) { DEL_PART; } // DEBUG TEST
 
769
 
 
770
#if LBM_INCLUDE_TESTSOLVERS==1
 
771
                        /*LbmFloat prob = 1.0;
 
772
                        // vanishing
 
773
                        prob = (rand()/(RAND_MAX+1.0));
 
774
                        // increse del prob. up to max height by given factor
 
775
                        const int fhCheckh = (int)(mpTest->mFluidHeight*1.25);
 
776
                        const LbmFloat fhProbh = 25.;
 
777
                        if((useff)&&(k>fhCheckh)) {
 
778
                                LbmFloat fac = (LbmFloat)(k-fhCheckh)/(LbmFloat)(fhProbh*(mLevel[mMaxRefine].lSizez-fhCheckh));
 
779
                                prob /= fac; 
 
780
                        }
 
781
                        //if(prob<mLevel[mMaxRefine].timestep*0.1) DEL_PART;
 
782
                        // forced vanishing
 
783
                        //? if(k>this->mSizez*3/4) {    if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
 
784
                        // */
 
785
 
 
786
#else // LBM_INCLUDE_TESTSOLVERS==1
 
787
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
788
 
 
789
                        if(p->getStatus()&PART_IN) { // IN
 
790
                                //if((k<cutval)||(k>this->mSizez-1-cutval)) DEL_PART; 
 
791
                                if(k<cutval) DEL_PART; 
 
792
                                if(k>this->mSizez-1-cutval) { P_CHANGETYPE(p, PART_DROP ); continue; } // create new drop
 
793
                                //ntlVec3Gfx v = getVelocityAt(i,j,k);
 
794
                                rho = vx = vy = vz = 0.0;
 
795
 
 
796
                                // define from particletracer.h
 
797
#if MOVE_FLOATS==1
 
798
                                const int DEPTH_AVG=3; // only average interface vels
 
799
                                int ccnt=0;
 
800
                                for(int kk=0;kk<DEPTH_AVG;kk+=1) {
 
801
                                //for(int kk=1;kk<DEPTH_AVG;kk+=1) {
 
802
                                        if((k-kk)<1) continue;
 
803
                                        //if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue;
 
804
                                        if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue;
 
805
                                        ccnt++;
 
806
                                        FORDF0{
 
807
                                                LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l);
 
808
                                                df[l] = cdf;
 
809
                                                //rho += cdf; 
 
810
                                                vx  += (this->dfDvecX[l]*cdf); 
 
811
                                                vy  += (this->dfDvecY[l]*cdf);  
 
812
                                                vz  += (this->dfDvecZ[l]*cdf);  
 
813
                                        }
 
814
                                }
 
815
                                if(ccnt) {
 
816
                                // use halved surface velocity (todo, use omega instead)
 
817
                                vx /=(LbmFloat)(ccnt * 2.0); // half xy speed! value2
 
818
                                vy /=(LbmFloat)(ccnt * 2.0);
 
819
                                vz /=(LbmFloat)(ccnt); }
 
820
#else // MOVE_FLOATS==1
 
821
                                vx=vy=0.; //p->setVel(ntlVec3Gfx(0.) ); // static_float
 
822
#endif // MOVE_FLOATS==1
 
823
                                vx += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
 
824
                                vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
 
825
 
 
826
                                //bool delfloat = false;
 
827
                                if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
 
828
                                        // in fluid cell
 
829
                                        /*if((1) && (k<this->mSizez-3) && 
 
830
                                                        (
 
831
                                                          TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFInter ) ||
 
832
                                                          TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) 
 
833
                                                        ) ) {
 
834
                                                vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
 
835
                                                if(vz<0.0) vz=0.0;
 
836
                                        } else delfloat=true;
 
837
                                        // keep below obstacles
 
838
                                        if((delfloat) && (TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFBnd )) ) {
 
839
                                                delfloat=false; vz=0.0;
 
840
                                        } // */
 
841
                                        vz = p->getVel()[2]-1.0*mLevel[mMaxRefine].gravity[2]; // simply rise...
 
842
                                        if(vz<0.) vz=0.;
 
843
                                } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFBnd ) ) {
 
844
                                        //vz = p->getVel()[2]-0.2; // force downwards movement, move below obstacle...
 
845
                                        vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall...
 
846
                                        if(vz>0.) vz=0.;
 
847
                                } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
 
848
                                        // keep in interface , one grid cell offset is added in part. gen
 
849
                                } else { //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { // shipt?, was all else before
 
850
                                        vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall...
 
851
                                        if(vz>0.) vz=0.;
 
852
                                        // check if above inter, remove otherwise
 
853
                                        /*if((1) && (k>2) && (
 
854
                                                                TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) ||
 
855
                                                                TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) 
 
856
                                                                ) ) {
 
857
                                                vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
 
858
                                                if(vz>0.0) vz=0.0;
 
859
                                        } else delfloat=true; // */
 
860
                                        //P_CHANGETYPE(p, PART_DROP ); continue; // create new drop
 
861
                                }
 
862
                                //if(delfloat) DEL_PART;
 
863
                                /*
 
864
                                // move down from empty
 
865
                                else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
 
866
                                        vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
 
867
                                        if(vz>0.0) vz=0.0;
 
868
                                        //DEL_PART; // ????
 
869
                                } else  {        DEL_PART; } // */
 
870
                                //vz = 0.0; // DEBUG
 
871
 
 
872
                                p->setVel( vec2G( ntlVec3Gfx(vx,vy,vz) ) ); //?
 
873
                                //p->setVel( vec2G(v)*0.75 + p->getVel()*0.25 ); //?
 
874
                                p->advanceVel();
 
875
                                //errMsg("PIT","IN pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
876
                        } else {
 
877
#if LBM_INCLUDE_TESTSOLVERS==1
 
878
                                if(useff) { mpTest->handleParticle(p, i,j,k); }
 
879
                                else DEL_PART; 
 
880
#else // LBM_INCLUDE_TESTSOLVERS==1
 
881
                                DEL_PART; 
 
882
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
883
                                //errMsg("PIT","OUT pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
 
884
                        }
 
885
                                
 
886
                        // additional bnd jitter
 
887
                        if((0) && (useff) && (p->getLifeTime()<3.*mLevel[mMaxRefine].timestep)) {
 
888
                                // use half butoff border 1/8
 
889
                                int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5);
 
890
                                if(maxdw<3) maxdw=3;
 
891
#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
 
892
#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) 
 
893
                                //if(ABS(i-(               cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
 
894
                                //if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5)   ,0.,0.); }
 
895
                                if((j>=0)&&(j<=this->mSizey-1)) {
 
896
                                        if(ABS(i-(               cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(               cutval))), 0.,0.); }
 
897
                                        if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); }
 
898
                                }
 
899
//#undef FLOAT_JITTER_BND
 
900
#undef FLOAT_JITTBNDRAND
 
901
                                //if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval)
 
902
                        }
 
903
                }  // PART_FLOAT
 
904
                
 
905
                // unknown particle type        
 
906
                else {
 
907
                        errMsg("LbmFsgrSolver::advanceParticles","PIT pit invalid type!? "<<p->getStatus() );
 
908
                }
367
909
  }
368
 
#endif // ELBEEM_PLUGIN
369
 
}
370
 
 
 
910
        myTime_t parttend = getTime(); 
 
911
        debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<", #particles:"<<mpParticles->getNumParticles() , 10 );
 
912
}
 
913
 
 
914
void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
 
915
        int workSet = mLevel[mMaxRefine].setCurr;
 
916
        std::ostringstream name;
 
917
 
 
918
        // debug - raw dump of ffrac values
 
919
        if(getenv("ELBEEM_RAWDEBUGDUMP")) {
 
920
                //name <<"fill_" << this->mStepCnt <<".dump";
 
921
                name << outfilename<< frameNrStr <<".dump";
 
922
                FILE *file = fopen(name.str().c_str(),"w");
 
923
                if(file) {
 
924
 
 
925
                        for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
 
926
                                for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
 
927
                                        for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
 
928
                                                float val = 0.;
 
929
                                                if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
 
930
                                                        val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
 
931
                                                        if(val<0.) val=0.;
 
932
                                                        if(val>1.) val=1.;
 
933
                                                }
 
934
                                                if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
 
935
                                                fprintf(file, "%f ",val); // text
 
936
                                                //errMsg("W", PRINT_IJK<<" val:"<<val);
 
937
                                        }
 
938
                                        fprintf(file, "\n"); // text
 
939
                                }
 
940
                                fprintf(file, "\n"); // text
 
941
                        }
 
942
                        fclose(file);
 
943
 
 
944
                } // file
 
945
        } // */
 
946
        if(getenv("ELBEEM_BINDEBUGDUMP")) {
 
947
                name << outfilename<< frameNrStr <<".bdump";
 
948
                FILE *file = fopen(name.str().c_str(),"w");
 
949
                if(file) {
 
950
                        for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
 
951
                                for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
 
952
                                        for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
 
953
                                                float val = 0.;
 
954
                                                if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
 
955
                                                        val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
 
956
                                                        if(val<0.) val=0.;
 
957
                                                        if(val>1.) val=1.;
 
958
                                                }
 
959
                                                if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
 
960
                                                fwrite( &val, sizeof(val), 1, file); // binary
 
961
                                        }
 
962
                                }
 
963
                        }
 
964
                        fclose(file);
 
965
                } // file
 
966
        }
 
967
 
 
968
        dumptype = 0; frameNr = 0; // get rid of warning
 
969
}
371
970
 
372
971
/*****************************************************************************/
373
972
/*! internal quick print function (for debugging) */
374
973
/*****************************************************************************/
375
 
template<class D>
376
974
void 
377
 
LbmFsgrSolver<D>::printLbmCell(int level, int i, int j, int k, int set) {
 
975
LbmFsgrSolver::printLbmCell(int level, int i, int j, int k, int set) {
378
976
        stdCellId *newcid = new stdCellId;
379
977
        newcid->level = level;
380
978
        newcid->x = i;
385
983
        debugPrintNodeInfo( newcid, set );
386
984
        delete newcid;
387
985
}
388
 
template<class D>
389
986
void 
390
 
LbmFsgrSolver<D>::debugMarkCellCall(int level, int vi,int vj,int vk) {
 
987
LbmFsgrSolver::debugMarkCellCall(int level, int vi,int vj,int vk) {
391
988
        stdCellId *newcid = new stdCellId;
392
989
        newcid->level = level;
393
990
        newcid->x = vi;
394
991
        newcid->y = vj;
395
992
        newcid->z = vk;
396
 
        addCellToMarkedList( newcid );
 
993
        this->addCellToMarkedList( newcid );
397
994
}
398
995
 
399
996
                
414
1011
#define CID_EY (int)( (mLevel[cid->level].lSizey-1) * guiRoiEY )
415
1012
#define CID_EZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiEZ )
416
1013
 
417
 
template<class D>
418
1014
CellIdentifierInterface* 
419
 
LbmFsgrSolver<D>::getFirstCell( ) {
 
1015
LbmFsgrSolver::getFirstCell( ) {
420
1016
        int level = mMaxRefine;
421
1017
 
422
1018
#if LBMDIM==3
434
1030
        return cid;
435
1031
}
436
1032
 
437
 
template<class D>
438
 
typename LbmFsgrSolver<D>::stdCellId* 
439
 
LbmFsgrSolver<D>::convertBaseCidToStdCid( CellIdentifierInterface* basecid) {
 
1033
LbmFsgrSolver::stdCellId* 
 
1034
LbmFsgrSolver::convertBaseCidToStdCid( CellIdentifierInterface* basecid) {
440
1035
        //stdCellId *cid = dynamic_cast<stdCellId*>( basecid );
441
1036
        stdCellId *cid = (stdCellId*)( basecid );
442
1037
        return cid;
443
1038
}
444
1039
 
445
 
template<class D>
446
 
void 
447
 
LbmFsgrSolver<D>::advanceCell( CellIdentifierInterface* basecid) {
 
1040
void LbmFsgrSolver::advanceCell( CellIdentifierInterface* basecid) {
448
1041
        stdCellId *cid = convertBaseCidToStdCid(basecid);
449
1042
        if(cid->getEnd()) return;
450
1043
 
467
1060
        //debugOut(" ADa "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10);
468
1061
}
469
1062
 
470
 
template<class D>
471
 
bool 
472
 
LbmFsgrSolver<D>::noEndCell( CellIdentifierInterface* basecid) {
 
1063
bool LbmFsgrSolver::noEndCell( CellIdentifierInterface* basecid) {
473
1064
        stdCellId *cid = convertBaseCidToStdCid(basecid);
474
1065
        return (!cid->getEnd());
475
1066
}
476
1067
 
477
 
template<class D>
478
 
void 
479
 
LbmFsgrSolver<D>::deleteCellIterator( CellIdentifierInterface** cid ) {
 
1068
void LbmFsgrSolver::deleteCellIterator( CellIdentifierInterface** cid ) {
480
1069
        delete *cid;
481
1070
        *cid = NULL;
482
1071
}
483
1072
 
484
 
template<class D>
485
 
CellIdentifierInterface* 
486
 
LbmFsgrSolver<D>::getCellAt( ntlVec3Gfx pos ) {
 
1073
CellIdentifierInterface* LbmFsgrSolver::getCellAt( ntlVec3Gfx pos ) {
487
1074
        //int cellok = false;
488
 
        pos -= (D::mvGeoStart);
 
1075
        pos -= (this->mvGeoStart);
489
1076
 
490
1077
        LbmFloat mmaxsize = mLevel[mMaxRefine].nodeSize;
491
1078
        for(int level=mMaxRefine; level>=0; level--) { // finest first
492
1079
        //for(int level=0; level<=mMaxRefine; level++) { // coarsest first
493
1080
                LbmFloat nsize = mLevel[level].nodeSize;
494
1081
                int x,y,z;
495
 
                //LbmFloat nsize = getCellSize(NULL)[0]*2.0;
496
 
                x = (int)((pos[0]-0.5*mmaxsize) / nsize );
497
 
                y = (int)((pos[1]-0.5*mmaxsize) / nsize );
498
 
                z = (int)((pos[2]-0.5*mmaxsize) / nsize );
499
 
                if(D::cDimension==2) z = 0;
 
1082
                // CHECK +- maxsize?
 
1083
                x = (int)((pos[0]+0.5*mmaxsize) / nsize );
 
1084
                y = (int)((pos[1]+0.5*mmaxsize) / nsize );
 
1085
                z = (int)((pos[2]+0.5*mmaxsize) / nsize );
 
1086
                if(LBMDIM==2) z = 0;
500
1087
 
501
1088
                // double check...
502
 
                //int level = mMaxRefine;
503
1089
                if(x<0) continue;
504
1090
                if(y<0) continue;
505
1091
                if(z<0) continue;
518
1104
                newcid->x = x;
519
1105
                newcid->y = y;
520
1106
                newcid->z = z;
521
 
                //errMsg("cellAt",D::mName<<" "<<pos<<" l"<<level<<":"<<x<<","<<y<<","<<z<<" "<<convertCellFlagType2String(RFLAG(level, x,y,z, mLevel[level].setCurr)) );
 
1107
                //errMsg("cellAt",this->mName<<" "<<pos<<" l"<<level<<":"<<x<<","<<y<<","<<z<<" "<<convertCellFlagType2String(RFLAG(level, x,y,z, mLevel[level].setCurr)) );
522
1108
                return newcid;
523
1109
        }
524
1110
 
528
1114
 
529
1115
// INFO functions
530
1116
 
531
 
template<class D>
532
 
int      
533
 
LbmFsgrSolver<D>::getCellSet      ( CellIdentifierInterface* basecid) {
 
1117
int      LbmFsgrSolver::getCellSet      ( CellIdentifierInterface* basecid) {
534
1118
        stdCellId *cid = convertBaseCidToStdCid(basecid);
535
1119
        return mLevel[cid->level].setCurr;
536
1120
        //return mLevel[cid->level].setOther;
537
1121
}
538
1122
 
539
 
template<class D>
540
 
int      
541
 
LbmFsgrSolver<D>::getCellLevel    ( CellIdentifierInterface* basecid) {
 
1123
int      LbmFsgrSolver::getCellLevel    ( CellIdentifierInterface* basecid) {
542
1124
        stdCellId *cid = convertBaseCidToStdCid(basecid);
543
1125
        return cid->level;
544
1126
}
545
1127
 
546
 
template<class D>
547
 
ntlVec3Gfx   
548
 
LbmFsgrSolver<D>::getCellOrigin   ( CellIdentifierInterface* basecid) {
 
1128
ntlVec3Gfx   LbmFsgrSolver::getCellOrigin   ( CellIdentifierInterface* basecid) {
549
1129
        ntlVec3Gfx ret;
550
1130
 
551
1131
        stdCellId *cid = convertBaseCidToStdCid(basecid);
552
1132
        ntlVec3Gfx cs( mLevel[cid->level].nodeSize );
553
 
        if(D::cDimension==2) { cs[2] = 0.0; }
 
1133
        if(LBMDIM==2) { cs[2] = 0.0; }
554
1134
 
555
 
        if(D::cDimension==2) {
556
 
                ret =(D::mvGeoStart -(cs*0.5) + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], (D::mvGeoEnd[2]-D::mvGeoStart[2])*0.5 )
 
1135
        if(LBMDIM==2) {
 
1136
                ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 )
557
1137
                                + ntlVec3Gfx(0.0,0.0,cs[1]*-0.25)*cid->level )
558
1138
                        +getCellSize(basecid);
559
1139
        } else {
560
 
                ret =(D::mvGeoStart -(cs*0.5) + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], cid->z *cs[2] ))
 
1140
                ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], cid->z *cs[2] ))
561
1141
                        +getCellSize(basecid);
562
1142
        }
563
1143
        return (ret);
564
1144
}
565
1145
 
566
 
template<class D>
567
 
ntlVec3Gfx   
568
 
LbmFsgrSolver<D>::getCellSize     ( CellIdentifierInterface* basecid) {
 
1146
ntlVec3Gfx   LbmFsgrSolver::getCellSize     ( CellIdentifierInterface* basecid) {
569
1147
        // return half size
570
1148
        stdCellId *cid = convertBaseCidToStdCid(basecid);
571
1149
        ntlVec3Gfx retvec( mLevel[cid->level].nodeSize * 0.5 );
572
1150
        // 2d display as rectangles
573
 
        if(D::cDimension==2) { retvec[2] = 0.0; }
 
1151
        if(LBMDIM==2) { retvec[2] = 0.0; }
574
1152
        return (retvec);
575
1153
}
576
1154
 
577
 
template<class D>
578
 
LbmFloat 
579
 
LbmFsgrSolver<D>::getCellDensity  ( CellIdentifierInterface* basecid,int set) {
 
1155
LbmFloat LbmFsgrSolver::getCellDensity  ( CellIdentifierInterface* basecid,int set) {
580
1156
        stdCellId *cid = convertBaseCidToStdCid(basecid);
581
1157
 
582
1158
        LbmFloat rho = 0.0;
583
 
        //FORDF0 { rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l); } // ORG
584
 
        //return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].stepsize) +1.0; // ORG
585
 
        if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) { // test
 
1159
        FORDF0 { rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l); } // ORG
 
1160
        return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep) +1.0; // ORG
 
1161
        /*if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) { // test
586
1162
                LbmFloat ux,uy,uz;
587
1163
                ux=uy=uz= 0.0;
588
1164
                int lev = cid->level;
589
1165
                LbmFloat df[27], feqOld[27];
590
1166
                FORDF0 {
591
1167
                        rho += QCELL(lev, cid->x,cid->y,cid->z, set, l);
592
 
                        ux += D::dfDvecX[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
593
 
                        uy += D::dfDvecY[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
594
 
                        uz += D::dfDvecZ[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
 
1168
                        ux += this->dfDvecX[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
 
1169
                        uy += this->dfDvecY[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
 
1170
                        uz += this->dfDvecZ[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
595
1171
                        df[l] = QCELL(lev, cid->x,cid->y,cid->z, set, l);
596
1172
                }
597
1173
                FORDF0 {
598
 
                        feqOld[l] = D::getCollideEq(l, rho,ux,uy,uz); 
 
1174
                        feqOld[l] = getCollideEq(l, rho,ux,uy,uz); 
599
1175
                }
600
 
                const LbmFloat Qo = D::getLesNoneqTensorCoeff(df,feqOld);
601
 
                //const LbmFloat modOmega = D::getLesOmega(mLevel[lev].omega, mLevel[lev].lcsmago,Qo);
 
1176
                // debugging mods
 
1177
                //const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feqOld);
 
1178
                //const LbmFloat modOmega = this->getLesOmega(mLevel[lev].omega, mLevel[lev].lcsmago,Qo);
602
1179
                //rho = (2.0-modOmega) *25.0;
603
 
                rho = Qo*100.0;
 
1180
                //rho = Qo*100.0;
604
1181
                //if(cid->x==24){ errMsg("MODOMT"," at "<<PRINT_VEC(cid->x,cid->y,cid->z)<<" = "<<rho<<" "<<Qo); }
605
1182
                //else{ rho=0.0; }
606
 
        } // test
607
 
        return rho; // test
 
1183
        } // test 
 
1184
        return rho; // test */
608
1185
}
609
1186
 
610
 
template<class D>
611
 
LbmVec   
612
 
LbmFsgrSolver<D>::getCellVelocity ( CellIdentifierInterface* basecid,int set) {
 
1187
LbmVec   LbmFsgrSolver::getCellVelocity ( CellIdentifierInterface* basecid,int set) {
613
1188
        stdCellId *cid = convertBaseCidToStdCid(basecid);
614
1189
 
 
1190
        // skip non-fluid cells
 
1191
        if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&(CFFluid|CFInter)) {
 
1192
                // ok go on...
 
1193
        } else {
 
1194
                return LbmVec(0.0);
 
1195
        }
 
1196
 
615
1197
        LbmFloat ux,uy,uz;
616
1198
        ux=uy=uz= 0.0;
617
1199
        FORDF0 {
618
 
                ux += D::dfDvecX[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
619
 
                uy += D::dfDvecY[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
620
 
                uz += D::dfDvecZ[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
 
1200
                ux += this->dfDvecX[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
 
1201
                uy += this->dfDvecY[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
 
1202
                uz += this->dfDvecZ[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
621
1203
        }
622
1204
        LbmVec vel(ux,uy,uz);
623
1205
        // TODO fix...
624
 
        return (vel * mLevel[cid->level].simCellSize / mLevel[cid->level].stepsize * D::mDebugVelScale); // normal
 
1206
        return (vel * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep * this->mDebugVelScale); // normal
625
1207
}
626
1208
 
627
 
template<class D>
628
 
LbmFloat   
629
 
LbmFsgrSolver<D>::getCellDf( CellIdentifierInterface* basecid,int set, int dir) {
 
1209
LbmFloat   LbmFsgrSolver::getCellDf( CellIdentifierInterface* basecid,int set, int dir) {
630
1210
        stdCellId *cid = convertBaseCidToStdCid(basecid);
631
1211
        return QCELL(cid->level, cid->x,cid->y,cid->z, set, dir);
632
1212
}
633
 
template<class D>
634
 
LbmFloat   
635
 
LbmFsgrSolver<D>::getCellMass( CellIdentifierInterface* basecid,int set) {
 
1213
LbmFloat   LbmFsgrSolver::getCellMass( CellIdentifierInterface* basecid,int set) {
636
1214
        stdCellId *cid = convertBaseCidToStdCid(basecid);
637
1215
        return QCELL(cid->level, cid->x,cid->y,cid->z, set, dMass);
638
1216
}
639
 
template<class D>
640
 
LbmFloat   
641
 
LbmFsgrSolver<D>::getCellFill( CellIdentifierInterface* basecid,int set) {
 
1217
LbmFloat   LbmFsgrSolver::getCellFill( CellIdentifierInterface* basecid,int set) {
642
1218
        stdCellId *cid = convertBaseCidToStdCid(basecid);
643
1219
        if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
644
1220
        if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFFluid) return 1.0;
645
1221
        return 0.0;
646
1222
        //return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
647
1223
}
648
 
template<class D>
649
 
CellFlagType 
650
 
LbmFsgrSolver<D>::getCellFlag( CellIdentifierInterface* basecid,int set) {
 
1224
CellFlagType LbmFsgrSolver::getCellFlag( CellIdentifierInterface* basecid,int set) {
651
1225
        stdCellId *cid = convertBaseCidToStdCid(basecid);
652
1226
        return RFLAG(cid->level, cid->x,cid->y,cid->z, set);
653
1227
}
654
1228
 
655
 
template<class D>
656
 
LbmFloat 
657
 
LbmFsgrSolver<D>::getEquilDf( int l ) {
658
 
        return D::dfEquil[l];
 
1229
LbmFloat LbmFsgrSolver::getEquilDf( int l ) {
 
1230
        return this->dfEquil[l];
659
1231
}
660
1232
 
661
 
template<class D>
662
 
int 
663
 
LbmFsgrSolver<D>::getDfNum( ) {
664
 
        return D::cDfNum;
 
1233
 
 
1234
ntlVec3Gfx LbmFsgrSolver::getVelocityAt   (float xp, float yp, float zp) {
 
1235
        ntlVec3Gfx avgvel(0.0);
 
1236
        LbmFloat   avgnum = 0.;
 
1237
 
 
1238
        // taken from getCellAt!
 
1239
        const int level = mMaxRefine;
 
1240
        const int workSet = mLevel[level].setCurr;
 
1241
        const LbmFloat nsize = mLevel[level].nodeSize;
 
1242
        const int x = (int)((-this->mvGeoStart[0]+xp-0.5*nsize) / nsize );
 
1243
        const int y = (int)((-this->mvGeoStart[1]+yp-0.5*nsize) / nsize );
 
1244
        int       z = (int)((-this->mvGeoStart[2]+zp-0.5*nsize) / nsize );
 
1245
        if(LBMDIM==2) z=0;
 
1246
        //errMsg("DUMPVEL","p"<<PRINT_VEC(xp,yp,zp)<<" at "<<PRINT_VEC(x,y,z)<<" max"<<PRINT_VEC(mLevel[level].lSizex,mLevel[level].lSizey,mLevel[level].lSizez) );
 
1247
 
 
1248
        // return fluid/if/border cells
 
1249
        // search neighborhood, do smoothing
 
1250
        FORDF0{ 
 
1251
                const int i = x+this->dfVecX[l];
 
1252
                const int j = y+this->dfVecY[l];
 
1253
                const int k = z+this->dfVecZ[l];
 
1254
 
 
1255
                if( (i<0) || (j<0) || (k<0) 
 
1256
                 || (i>=mLevel[level].lSizex) 
 
1257
                 || (j>=mLevel[level].lSizey) 
 
1258
                 || (k>=mLevel[level].lSizez) ) continue;
 
1259
 
 
1260
                if( (RFLAG(level, i,j,k, mLevel[level].setCurr)&(CFFluid|CFInter)) ) {
 
1261
                        ntlVec3Gfx vel(0.0);
 
1262
                        LbmFloat *ccel = RACPNT(level, i,j,k ,workSet); // omp
 
1263
                        for(int n=1; n<this->cDfNum; n++) {
 
1264
                                vel[0]  += (this->dfDvecX[n]*RAC(ccel,n));
 
1265
                                vel[1]  += (this->dfDvecY[n]*RAC(ccel,n)); 
 
1266
                                vel[2]  += (this->dfDvecZ[n]*RAC(ccel,n)); 
 
1267
                        } 
 
1268
 
 
1269
                        avgvel += vel;
 
1270
                        avgnum += 1.0;
 
1271
                        if(l==0) { // center slightly more weight
 
1272
                                avgvel += vel; avgnum += 1.0;
 
1273
                        }
 
1274
                } // */
 
1275
        }
 
1276
 
 
1277
        if(avgnum>0.) {
 
1278
                ntlVec3Gfx retv = avgvel / avgnum;
 
1279
                retv *= nsize/mLevel[level].timestep;
 
1280
                // scale for current animation settings (frame time)
 
1281
                retv *= mpParam->getCurrentAniFrameTime();
 
1282
                //errMsg("DUMPVEL","t"<<mSimulationTime<<" at "<<PRINT_VEC(xp,yp,zp)<<" ret:"<<retv<<", avgv:"<<avgvel<<" n"<<avgnum<<" nsize"<<nsize<<" ts"<<mLevel[level].timestep<<" fr"<<mpParam->getCurrentAniFrameTime() );
 
1283
                return retv;
 
1284
        }
 
1285
        // no cells here...?
 
1286
        //errMsg("DUMPVEL"," at "<<PRINT_VEC(xp,yp,zp)<<" v"<<avgvel<<" n"<<avgnum<<" no vel !?");
 
1287
        return ntlVec3Gfx(0.);
665
1288
}
666
1289
 
667
1290
#if LBM_USE_GUI==1
668
1291
//! show simulation info (implement SimulationObject pure virtual func)
669
 
template<class D>
670
1292
void 
671
 
LbmFsgrSolver<D>::debugDisplay(fluidDispSettings *set){ 
672
 
        //lbmDebugDisplay< LbmFsgrSolver<D> >( set, this ); 
 
1293
LbmFsgrSolver::debugDisplay(int set){ 
 
1294
        //lbmDebugDisplay< LbmFsgrSolver >( set, this ); 
673
1295
        lbmDebugDisplay( set ); 
674
1296
}
675
1297
#endif
680
1302
#if FSGR_STRICT_DEBUG==1
681
1303
#define STRICT_EXIT *((int *)0)=0;
682
1304
 
683
 
template<class D>
684
 
int LbmFsgrSolver<D>::debLBMGI(int level, int ii,int ij,int ik, int is) {
 
1305
int LbmFsgrSolver::debLBMGI(int level, int ii,int ij,int ik, int is) {
685
1306
        if(level <  0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
686
1307
        if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
687
1308
 
696
1317
        return _LBMGI(level, ii,ij,ik, is);
697
1318
};
698
1319
 
699
 
template<class D>
700
 
CellFlagType& LbmFsgrSolver<D>::debRFLAG(int level, int xx,int yy,int zz,int set){
 
1320
CellFlagType& LbmFsgrSolver::debRFLAG(int level, int xx,int yy,int zz,int set){
701
1321
        return _RFLAG(level, xx,yy,zz,set);   
702
1322
};
703
1323
 
704
 
template<class D>
705
 
CellFlagType& LbmFsgrSolver<D>::debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) {
 
1324
CellFlagType& LbmFsgrSolver::debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) {
706
1325
        if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
707
1326
        // warning might access all spatial nbs
708
 
        if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
 
1327
        if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
709
1328
        return _RFLAG_NB(level, xx,yy,zz,set, dir);
710
1329
};
711
1330
 
712
 
template<class D>
713
 
CellFlagType& LbmFsgrSolver<D>::debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) {
 
1331
CellFlagType& LbmFsgrSolver::debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) {
714
1332
        if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
715
 
        if(dir>D::cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
 
1333
        if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
716
1334
        return _RFLAG_NBINV(level, xx,yy,zz,set, dir);
717
1335
};
718
1336
 
719
 
template<class D>
720
 
int LbmFsgrSolver<D>::debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
 
1337
int LbmFsgrSolver::debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
721
1338
        if(level <  0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
722
1339
        if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
723
1340
 
730
1347
        if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
731
1348
        if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
732
1349
        if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
733
 
        if(l>D::cDfNum){  // dFfrac is an exception
 
1350
        if(l>this->cDfNum){  // dFfrac is an exception
734
1351
                if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
735
1352
#if COMPRESSGRIDS==1
736
 
        //if((!D::mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug
 
1353
        //if((!this->mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug
737
1354
#endif // COMPRESSGRIDS==1
738
1355
        return _LBMQI(level, ii,ij,ik, is, l);
739
1356
};
740
1357
 
741
 
template<class D>
742
 
LbmFloat& LbmFsgrSolver<D>::debQCELL(int level, int xx,int yy,int zz,int set,int l) {
 
1358
LbmFloat& LbmFsgrSolver::debQCELL(int level, int xx,int yy,int zz,int set,int l) {
743
1359
        //errMsg("LbmStrict","debQCELL debug: l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" l"<<l<<" index"<<LBMGI(level, xx,yy,zz,set)); 
744
1360
        return _QCELL(level, xx,yy,zz,set,l);
745
1361
};
746
1362
 
747
 
template<class D>
748
 
LbmFloat& LbmFsgrSolver<D>::debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) {
 
1363
LbmFloat& LbmFsgrSolver::debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) {
749
1364
        if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
750
 
        if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
 
1365
        if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
751
1366
        return _QCELL_NB(level, xx,yy,zz,set, dir,l);
752
1367
};
753
1368
 
754
 
template<class D>
755
 
LbmFloat& LbmFsgrSolver<D>::debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) {
 
1369
LbmFloat& LbmFsgrSolver::debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) {
756
1370
        if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
757
 
        if(dir>D::cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
 
1371
        if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
758
1372
        return _QCELL_NBINV(level, xx,yy,zz,set, dir,l);
759
1373
};
760
1374
 
761
 
template<class D>
762
 
LbmFloat* LbmFsgrSolver<D>::debRACPNT(int level,  int ii,int ij,int ik, int is ) {
 
1375
LbmFloat* LbmFsgrSolver::debRACPNT(int level,  int ii,int ij,int ik, int is ) {
763
1376
        return _RACPNT(level, ii,ij,ik, is );
764
1377
};
765
1378
 
766
 
template<class D>
767
 
LbmFloat& LbmFsgrSolver<D>::debRAC(LbmFloat* s,int l) {
 
1379
LbmFloat& LbmFsgrSolver::debRAC(LbmFloat* s,int l) {
768
1380
        if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
769
1381
        if(l>dTotalNum){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } 
770
 
        //if(l>D::cDfNum){ // dFfrac is an exception 
 
1382
        //if(l>this->cDfNum){ // dFfrac is an exception 
771
1383
        //if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
772
1384
        return _RAC(s,l);
773
1385
};
785
1397
#include "../gui/gui_utilities.h"
786
1398
 
787
1399
//! display a single node
788
 
template<class D>
789
 
void LbmFsgrSolver<D>::debugDisplayNode(fluidDispSettings *dispset, CellIdentifierInterface* cell ) {
 
1400
void LbmFsgrSolver::debugDisplayNode(int dispset, CellIdentifierInterface* cell ) {
790
1401
        //debugOut(" DD: "<<cell->getAsString() , 10);
791
1402
        ntlVec3Gfx org      = this->getCellOrigin( cell );
792
1403
        ntlVec3Gfx halfsize = this->getCellSize( cell );
796
1407
        bool     showcell = true;
797
1408
        int      linewidth = 1;
798
1409
        ntlColor col(0.5);
799
 
        LbmFloat cscale = dispset->scale;
 
1410
        LbmFloat cscale = 1.0; //dispset->scale;
800
1411
 
801
1412
#define DRAWDISPCUBE(col,scale) \
802
1413
        {       glLineWidth( linewidth ); \
821
1432
        else
822
1433
        if(flag& CFFluid    )    { if(!guiShowFluid    ) return; }
823
1434
 
824
 
        switch(dispset->type) {
 
1435
        switch(dispset) {
825
1436
                case FLUIDDISPNothing: {
826
1437
                                showcell = false;
827
1438
                        } break;
947
1558
 
948
1559
//! debug display function
949
1560
//  D has to implement the CellIterator interface
950
 
template<class D>
951
 
void LbmFsgrSolver<D>::lbmDebugDisplay(fluidDispSettings *dispset) {
952
 
        //je nach solver...?
953
 
        if(!dispset->on) return;
 
1561
void LbmFsgrSolver::lbmDebugDisplay(int dispset) {
 
1562
        // DEBUG always display testdata
 
1563
#if LBM_INCLUDE_TESTSOLVERS==1
 
1564
        if(mUseTestdata){ mpTest->testDebugDisplay(dispset); }
 
1565
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
1566
        if(dispset<=FLUIDDISPNothing) return;
 
1567
        //if(!dispset->on) return;
954
1568
        glDisable( GL_LIGHTING ); // dont light lines
955
1569
 
956
 
        typename D::CellIdentifier cid = this->getFirstCell();
 
1570
#if LBM_INCLUDE_TESTSOLVERS==1
 
1571
        if((!mUseTestdata)|| (mUseTestdata)&&(mpTest->mDebugvalue1<=0.0)) {
 
1572
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
1573
 
 
1574
        LbmFsgrSolver::CellIdentifier cid = this->getFirstCell();
957
1575
        for(; this->noEndCell( cid );
958
1576
              this->advanceCell( cid ) ) {
959
1577
                this->debugDisplayNode(dispset, cid );
960
1578
        }
961
1579
        delete cid;
962
1580
 
 
1581
#if LBM_INCLUDE_TESTSOLVERS==1
 
1582
        } // 3d check
 
1583
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
1584
 
963
1585
        glEnable( GL_LIGHTING ); // dont light lines
964
1586
}
965
1587
 
966
1588
//! debug display function
967
1589
//  D has to implement the CellIterator interface
968
 
template<class D>
969
 
void LbmFsgrSolver<D>::lbmMarkedCellDisplay() {
970
 
        fluidDispSettings dispset;
 
1590
void LbmFsgrSolver::lbmMarkedCellDisplay() {
 
1591
        //fluidDispSettings dispset;
971
1592
        // trick - display marked cells as grid displa -> white, big
972
 
        dispset.type = FLUIDDISPGrid;
973
 
        dispset.on = true;
 
1593
        int dispset = FLUIDDISPGrid;
974
1594
        glDisable( GL_LIGHTING ); // dont light lines
975
1595
        
976
 
        typename D::CellIdentifier cid = this->markedGetFirstCell();
 
1596
        LbmFsgrSolver::CellIdentifier cid = this->markedGetFirstCell();
977
1597
        while(cid) {
978
 
                this->debugDisplayNode(&dispset, cid );
 
1598
                this->debugDisplayNode(dispset, cid );
979
1599
                cid = this->markedAdvanceCell();
980
1600
        }
981
1601
        delete cid;
986
1606
#endif // LBM_USE_GUI==1
987
1607
 
988
1608
//! display a single node
989
 
template<class D>
990
 
void LbmFsgrSolver<D>::debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet) {
 
1609
void LbmFsgrSolver::debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet) {
991
1610
                //string printInfo,
992
1611
                // force printing of one set? default = -1 = off
993
1612
  bool printDF     = false;
1037
1656
                debMsgStd("                  ",DM_MSG, "Printing set:"<<workset<<" orgSet:"<<set, 1);
1038
1657
                
1039
1658
                if(printDF) {
1040
 
                        for(int l=0; l<this->getDfNum(); l++) { // FIXME ??
 
1659
                        for(int l=0; l<LBM_DFNUM; l++) { // FIXME ??
1041
1660
                                debMsgStd("                  ",DM_MSG, "  Df"<<l<<": "<<this->getCellDf(cell,workset,l), 1);
1042
1661
                        }
1043
1662
                }
1054
1673
                if(printMass) {
1055
1674
                        debMsgStd("                  ",DM_MSG, "  Mss: "<<this->getCellMass(cell,workset), 1);
1056
1675
                }
 
1676
                // dirty... TODO fixme
 
1677
                debMsgStd("                  ",DM_MSG, "  Flx: "<<this->getCellDf(cell,workset,dFlux), 1);
1057
1678
        }
1058
1679
}
1059
1680
 
1060
 
#endif // !defined(__APPLE_CC__) || defined(LBM_FORCEINCLUDE)
1061
 
 
1062
 
/******************************************************************************
1063
 
 * instantiation
1064
 
 *****************************************************************************/
1065
 
#if ((!defined(__APPLE_CC__)) && (!defined(__INTEL_COMPILER))) && (!defined(LBM_FORCEINCLUDE))
1066
 
 
1067
 
#if LBMDIM==2
1068
 
#define LBM_INSTANTIATE LbmBGK2D
1069
 
#endif // LBMDIM==2
1070
 
#if LBMDIM==3
1071
 
#define LBM_INSTANTIATE LbmBGK3D
1072
 
#endif // LBMDIM==3
1073
 
 
1074
 
template class LbmFsgrSolver< LBM_INSTANTIATE >;
1075
 
 
1076
 
#endif // __APPLE_CC__ __INTEL_COMPILER
1077
1681