← 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_class.h

  • 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_class.h
15
15
 
16
16
#include "utilities.h"
17
17
#include "solver_interface.h"
18
 
#include "ntl_scene.h"
 
18
#include "ntl_ray.h"
19
19
#include <stdio.h>
20
20
 
21
21
#if PARALLEL==1
25
25
#define PARALLEL 0
26
26
#endif // PARALLEL
27
27
 
28
 
#ifndef LBMMODEL_DEFINED
29
 
// force compiler error!
30
 
ERROR - define model first!
31
 
#endif // LBMMODEL_DEFINED
32
 
 
33
28
 
34
29
// general solver setting defines
35
30
 
36
 
// default to 3dim
37
 
#ifndef LBMDIM
38
 
#define LBMDIM 3
39
 
#endif // LBMDIM
40
 
 
41
 
 
42
31
//! debug coordinate accesses and the like? (much slower)
 
32
//  might be enabled by compilation
 
33
#ifndef FSGR_STRICT_DEBUG
43
34
#define FSGR_STRICT_DEBUG 0
 
35
#endif // FSGR_STRICT_DEBUG
44
36
 
45
37
//! debug coordinate accesses and the like? (much slower)
46
38
#define FSGR_OMEGA_DEBUG 0
81
73
#define FSGR_MAXNOOFLEVELS 5
82
74
 
83
75
// enable/disable fine grid compression for finest level
 
76
// make sure this is same as useGridComp in calculateMemreqEstimate
84
77
#if LBMDIM==3
85
78
#define COMPRESSGRIDS 1
86
79
#else 
113
106
#endif
114
107
#endif
115
108
 
116
 
#if ELBEEM_BLENDER!=1
117
 
#include "solver_test.h"
118
 
#endif // ELBEEM_BLENDER==1
119
 
 
120
109
/*****************************************************************************/
121
110
/*! cell access classes */
122
 
template<typename D>
123
111
class UniformFsgrCellIdentifier : 
124
 
        public CellIdentifierInterface 
 
112
        public CellIdentifierInterface , public LbmCellContents
125
113
{
126
114
        public:
127
115
                //! which grid level?
137
125
                virtual string getAsString() {
138
126
                        std::ostringstream ret;
139
127
                        ret <<"{ i"<<x<<",j"<<y;
140
 
                        if(D::cDimension>2) ret<<",k"<<z;
 
128
                        if(LBMDIM>2) ret<<",k"<<z;
141
129
                        ret <<" }";
142
130
                        return ret.str();
143
131
                }
144
132
 
145
133
                virtual bool equal(CellIdentifierInterface* other) {
146
 
                        //UniformFsgrCellIdentifier<D> *cid = dynamic_cast<UniformFsgrCellIdentifier<D> *>( other );
147
 
                        UniformFsgrCellIdentifier<D> *cid = (UniformFsgrCellIdentifier<D> *)( other );
 
134
                        UniformFsgrCellIdentifier *cid = (UniformFsgrCellIdentifier *)( other );
148
135
                        if(!cid) return false;
149
136
                        if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true;
150
137
                        return false;
165
152
        //! size this level was advanced to 
166
153
        LbmFloat time;
167
154
        //! size of a single lbm step in time units on this level 
168
 
        LbmFloat stepsize;
 
155
        LbmFloat timestep;
169
156
        //! step count
170
157
        int lsteps;
171
158
        //! gravity force for this level
203
190
 
204
191
/*****************************************************************************/
205
192
/*! class for solving a LBM problem */
206
 
template<class D>
207
193
class LbmFsgrSolver : 
208
 
        public D // this means, the solver is a lbmData object and implements the lbmInterface
 
194
        public LbmSolverInterface // this means, the solver is a lbmData object and implements the lbmInterface
209
195
{
210
196
 
211
197
        public:
213
199
                LbmFsgrSolver();
214
200
                //! Destructor 
215
201
                virtual ~LbmFsgrSolver();
216
 
                //! id string of solver
217
 
                virtual string getIdString();
218
 
                //! dimension of solver
219
 
                virtual int getDimension();
220
202
 
221
203
                //! initilize variables fom attribute list 
222
204
                virtual void parseAttrList();
223
205
                //! Initialize omegas and forces on all levels (for init/timestep change)
224
206
                void initLevelOmegas();
225
 
                //! finish the init with config file values (allocate arrays...) 
226
 
                virtual bool initializeSolver(); //( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ );
 
207
 
 
208
                // multi step solver init
 
209
                /*! finish the init with config file values (allocate arrays...) */
 
210
                virtual bool initializeSolverMemory();
 
211
                /*! init solver arrays */
 
212
                virtual bool initializeSolverGrids();
 
213
                /*! prepare actual simulation start, setup viz etc */
 
214
                virtual bool initializeSolverPostinit();
 
215
 
 
216
                //! notify object that dump is in progress (e.g. for field dump) 
 
217
                virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
227
218
 
228
219
#if LBM_USE_GUI==1
229
220
                //! show simulation info (implement LbmSolverInterface pure virtual func)
230
 
                virtual void debugDisplay(fluidDispSettings *set);
 
221
                virtual void debugDisplay(int set);
231
222
#endif
232
223
                
233
 
                
234
224
                // implement CellIterator<UniformFsgrCellIdentifier> interface
235
 
                typedef UniformFsgrCellIdentifier<typename D::LbmCellContents> stdCellId;
 
225
                typedef UniformFsgrCellIdentifier stdCellId;
236
226
                virtual CellIdentifierInterface* getFirstCell( );
237
227
                virtual void advanceCell( CellIdentifierInterface* );
238
228
                virtual bool noEndCell( CellIdentifierInterface* );
249
239
                virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set);
250
240
                virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set);
251
241
                virtual LbmFloat   getEquilDf      ( int );
252
 
                virtual int        getDfNum        ( );
 
242
                virtual ntlVec3Gfx getVelocityAt   (float x, float y, float z);
253
243
                // convert pointers
254
244
                stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid);
255
245
 
264
254
                LBM_INLINED void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
265
255
                LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
266
256
                LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
 
257
                //! interpolate velocity and density at a given position
 
258
                void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz);
267
259
 
268
260
                /*! perform a single LBM step */
269
261
                void stepMain();
278
270
                /* simulation object interface, just calls stepMain */
279
271
                virtual void step();
280
272
                /*! init particle positions */
281
 
                virtual int initParticles(ParticleTracer *partt);
 
273
                virtual int initParticles();
282
274
                /*! move all particles */
283
 
                virtual void advanceParticles(ParticleTracer *partt );
 
275
                virtual void advanceParticles();
284
276
 
285
277
 
286
278
                /*! debug object display (used e.g. for preview surface) */
288
280
        
289
281
                // gui/output debugging functions
290
282
#if LBM_USE_GUI==1
291
 
                virtual void debugDisplayNode(fluidDispSettings *dispset, CellIdentifierInterface* cell );
292
 
                virtual void lbmDebugDisplay(fluidDispSettings *dispset);
 
283
                virtual void debugDisplayNode(int dispset, CellIdentifierInterface* cell );
 
284
                virtual void lbmDebugDisplay(int dispset);
293
285
                virtual void lbmMarkedCellDisplay();
294
286
#endif // LBM_USE_GUI==1
295
287
                virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1);
298
290
                void prepareVisualization( void );
299
291
 
300
292
                /*! type for cells */
301
 
                typedef typename D::LbmCell LbmCell;
 
293
                //typedef typename this->LbmCell LbmCell;
302
294
                
303
295
        protected:
304
296
 
311
303
                void adaptTimestep();
312
304
                //! init mObjectSpeeds for current parametrization
313
305
                void recalculateObjectSpeeds();
 
306
                //! init moving obstacles for next sim step sim 
 
307
                void initMovingObstacles(bool staticInit);
314
308
                //! flag reinit step - always works on finest grid!
315
309
                void reinitFlags( int workSet );
316
310
                //! mass dist weights
347
341
                
348
342
                //! use time adaptivity? 
349
343
                bool mTimeAdap;
350
 
                //! domain boundary free/no slip type
351
 
                string mDomainBound;
352
 
                //! part slip value for domain
353
 
                LbmFloat mDomainPartSlipValue;
 
344
                //! force smaller timestep for next LBM step? (eg for mov obj)
 
345
                bool mForceTimeStepReduce;
354
346
 
355
347
                //! fluid vol height
356
348
                LbmFloat mFVHeight;
364
356
                //! smoother surface initialization?
365
357
                int mInitSurfaceSmoothing;
366
358
 
 
359
                //! lock time step down switching
367
360
                int mTimestepReduceLock;
 
361
                //! count no. of switches
368
362
                int mTimeSwitchCounts;
 
363
                // only switch of maxvel is higher for several steps...
 
364
                int mTimeMaxvelStepCnt;
 
365
 
369
366
                //! total simulation time so far 
370
 
                LbmFloat mSimulationTime;
 
367
                LbmFloat mSimulationTime, mLastSimTime;
371
368
                //! smallest and largest step size so far 
372
 
                LbmFloat mMinStepTime, mMaxStepTime;
 
369
                LbmFloat mMinTimestep, mMaxTimestep;
373
370
                //! track max. velocity
374
371
                LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen;
375
372
 
403
400
                vector<LbmVec> mObjectSpeeds;
404
401
                //! partslip bc. values for obstacle boundary conditions
405
402
                vector<LbmFloat> mObjectPartslips;
 
403
                //! moving object mass boundary condition values
 
404
                vector<LbmFloat> mObjectMassMovnd;
 
405
 
 
406
                //! permanent movobj vert storage
 
407
          vector<ntlVec3Gfx>  mMOIVertices;
 
408
        vector<ntlVec3Gfx>  mMOIVerticesOld;
406
409
 
407
410
                //! get isofield weights
408
411
                int mIsoWeightMethod;
442
445
                int mDisableStandingFluidInit;
443
446
                //! debug function to force tadap syncing
444
447
                int mForceTadapRefine;
445
 
 
446
 
#ifndef ELBEEM_BLENDER
447
 
                // test functions
448
 
                bool mUseTestdata;
449
 
                LbmTestdata *mpTest;
450
 
                void initTestdata();
451
 
                void destroyTestdata();
452
 
                void handleTestdata();
453
 
                void exportTestdata();
454
 
                ParticleTracer *mpParticles;
455
 
#endif // ELBEEM_BLENDER==1
 
448
                //! border cutoff value
 
449
                int mCutoff;
456
450
 
457
451
                // strict debug interface
458
452
#               if FSGR_STRICT_DEBUG==1
467
461
                LbmFloat* debRACPNT(int level,  int ii,int ij,int ik, int is );
468
462
                LbmFloat& debRAC(LbmFloat* s,int l);
469
463
#               endif // FSGR_STRICT_DEBUG==1
 
464
 
 
465
                bool mUseTestdata;
 
466
 
 
467
        public: // former LbmModelLBGK  functions
 
468
                // relaxation funtions - implemented together with relax macros
 
469
                static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz);
 
470
                static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz);
 
471
                inline LbmFloat getLesNoneqTensorCoeff( LbmFloat df[],                          LbmFloat feq[] );
 
472
                inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo);
 
473
                inline void collideArrays( int i, int j, int k, // position - more for debugging
 
474
                                LbmFloat df[], LbmFloat &outrho, // out only!
 
475
                                // velocity modifiers (returns actual velocity!)
 
476
                                LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
 
477
                                LbmFloat omega, LbmVec gravity, LbmFloat csmago, 
 
478
                                LbmFloat *newOmegaRet, LbmFloat *newQoRet);
 
479
 
 
480
 
 
481
                // former LBM models
 
482
        public:
 
483
 
 
484
//! shorten static const definitions
 
485
#define STCON static const
 
486
 
 
487
#if LBMDIM==3
 
488
                
 
489
                //! id string of solver
 
490
                virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D3Q19]"); }
 
491
 
 
492
                //! how many dimensions? UNUSED? replace by LBMDIM?
 
493
                STCON int cDimension;
 
494
 
 
495
                // Wi factors for collide step 
 
496
                STCON LbmFloat cCollenZero;
 
497
                STCON LbmFloat cCollenOne;
 
498
                STCON LbmFloat cCollenSqrtTwo;
 
499
 
 
500
                //! threshold value for filled/emptied cells 
 
501
                STCON LbmFloat cMagicNr2;
 
502
                STCON LbmFloat cMagicNr2Neg;
 
503
                STCON LbmFloat cMagicNr;
 
504
                STCON LbmFloat cMagicNrNeg;
 
505
 
 
506
                //! size of a single set of distribution functions 
 
507
                STCON int    cDfNum;
 
508
                //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
 
509
                STCON int    cDirNum;
 
510
 
 
511
                //! distribution functions directions 
 
512
                typedef enum {
 
513
                         cDirInv=  -1,
 
514
                         cDirC  =  0,
 
515
                         cDirN  =  1,
 
516
                         cDirS  =  2,
 
517
                         cDirE  =  3,
 
518
                         cDirW  =  4,
 
519
                         cDirT  =  5,
 
520
                         cDirB  =  6,
 
521
                         cDirNE =  7,
 
522
                         cDirNW =  8,
 
523
                         cDirSE =  9,
 
524
                         cDirSW = 10,
 
525
                         cDirNT = 11,
 
526
                         cDirNB = 12,
 
527
                         cDirST = 13,
 
528
                         cDirSB = 14,
 
529
                         cDirET = 15,
 
530
                         cDirEB = 16,
 
531
                         cDirWT = 17,
 
532
                         cDirWB = 18
 
533
                } dfDir;
 
534
 
 
535
                /* Vector Order 3D:
 
536
                 *  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
 
537
                 *  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
 
538
                 *  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
 
539
                 *  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1
 
540
                 */
 
541
 
 
542
                /*! name of the dist. function 
 
543
                         only for nicer output */
 
544
                STCON char* dfString[ 19 ];
 
545
 
 
546
                /*! index of normal dist func, not used so far?... */
 
547
                STCON int dfNorm[ 19 ];
 
548
 
 
549
                /*! index of inverse dist func, not fast, but useful... */
 
550
                STCON int dfInv[ 19 ];
 
551
 
 
552
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
553
                STCON int dfRefX[ 19 ];
 
554
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
555
                STCON int dfRefY[ 19 ];
 
556
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
557
                STCON int dfRefZ[ 19 ];
 
558
 
 
559
                /*! dist func vectors */
 
560
                STCON int dfVecX[ 27 ];
 
561
                STCON int dfVecY[ 27 ];
 
562
                STCON int dfVecZ[ 27 ];
 
563
 
 
564
                /*! arrays as before with doubles */
 
565
                STCON LbmFloat dfDvecX[ 27 ];
 
566
                STCON LbmFloat dfDvecY[ 27 ];
 
567
                STCON LbmFloat dfDvecZ[ 27 ];
 
568
 
 
569
                /*! principal directions */
 
570
                STCON int princDirX[ 2*3 ];
 
571
                STCON int princDirY[ 2*3 ];
 
572
                STCON int princDirZ[ 2*3 ];
 
573
 
 
574
                /*! vector lengths */
 
575
                STCON LbmFloat dfLength[ 19 ];
 
576
 
 
577
                /*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
 
578
                static LbmFloat dfEquil[ 19 ];
 
579
 
 
580
                /*! arrays for les model coefficients */
 
581
                static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
 
582
                static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
 
583
 
 
584
#else // end LBMDIM==3 , LBMDIM==2
 
585
                
 
586
                //! id string of solver
 
587
                virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D2Q9]"); }
 
588
 
 
589
                //! how many dimensions?
 
590
                STCON int cDimension;
 
591
 
 
592
                //! Wi factors for collide step 
 
593
                STCON LbmFloat cCollenZero;
 
594
                STCON LbmFloat cCollenOne;
 
595
                STCON LbmFloat cCollenSqrtTwo;
 
596
 
 
597
                //! threshold value for filled/emptied cells 
 
598
                STCON LbmFloat cMagicNr2;
 
599
                STCON LbmFloat cMagicNr2Neg;
 
600
                STCON LbmFloat cMagicNr;
 
601
                STCON LbmFloat cMagicNrNeg;
 
602
 
 
603
                //! size of a single set of distribution functions 
 
604
                STCON int    cDfNum;
 
605
                STCON int    cDirNum;
 
606
 
 
607
                //! distribution functions directions 
 
608
                typedef enum {
 
609
                         cDirInv=  -1,
 
610
                         cDirC  =  0,
 
611
                         cDirN  =  1,
 
612
                         cDirS  =  2,
 
613
                         cDirE  =  3,
 
614
                         cDirW  =  4,
 
615
                         cDirNE =  5,
 
616
                         cDirNW =  6,
 
617
                         cDirSE =  7,
 
618
                         cDirSW =  8
 
619
                } dfDir;
 
620
 
 
621
                /* Vector Order 2D:
 
622
                 * 0  1 2  3  4  5  6 7  8
 
623
                 * 0, 0,0, 1,-1, 1,-1,1,-1 
 
624
                 * 0, 1,-1, 0,0, 1,1,-1,-1  */
 
625
 
 
626
                /* name of the dist. function 
 
627
                         only for nicer output */
 
628
                STCON char* dfString[ 9 ];
 
629
 
 
630
                /* index of normal dist func, not used so far?... */
 
631
                STCON int dfNorm[ 9 ];
 
632
 
 
633
                /* index of inverse dist func, not fast, but useful... */
 
634
                STCON int dfInv[ 9 ];
 
635
 
 
636
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
637
                STCON int dfRefX[ 9 ];
 
638
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
639
                STCON int dfRefY[ 9 ];
 
640
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
641
                STCON int dfRefZ[ 9 ];
 
642
 
 
643
                /* dist func vectors */
 
644
                STCON int dfVecX[ 9 ];
 
645
                STCON int dfVecY[ 9 ];
 
646
                /* Z, 2D values are all 0! */
 
647
                STCON int dfVecZ[ 9 ];
 
648
 
 
649
                /* arrays as before with doubles */
 
650
                STCON LbmFloat dfDvecX[ 9 ];
 
651
                STCON LbmFloat dfDvecY[ 9 ];
 
652
                /* Z, 2D values are all 0! */
 
653
                STCON LbmFloat dfDvecZ[ 9 ];
 
654
 
 
655
                /*! principal directions */
 
656
                STCON int princDirX[ 2*2 ];
 
657
                STCON int princDirY[ 2*2 ];
 
658
                STCON int princDirZ[ 2*2 ];
 
659
 
 
660
                /* vector lengths */
 
661
                STCON LbmFloat dfLength[ 9 ];
 
662
 
 
663
                /* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
 
664
                static LbmFloat dfEquil[ 9 ];
 
665
 
 
666
                /*! arrays for les model coefficients */
 
667
                static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
 
668
                static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
 
669
 
 
670
#endif  // LBMDIM==2
470
671
};
471
672
 
 
673
#undef STCON
472
674
 
473
675
 
474
676
/*****************************************************************************/
479
681
// cell mark debugging
480
682
#if FSGR_STRICT_DEBUG==10
481
683
#define debugMarkCell(lev,x,y,z) \
482
 
        errMsg("debugMarkCell",D::mName<<" step: "<<D::mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
 
684
        errMsg("debugMarkCell",this->mName<<" step: "<<this->mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
483
685
        debugMarkCellCall((lev),(x),(y),(z));
484
686
#else // FSGR_STRICT_DEBUG==1
485
687
#define debugMarkCell(lev,x,y,z) \
494
696
 
495
697
//! flag array acces macro
496
698
#define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ]
497
 
#define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set) ]
498
 
#define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set) ]
499
 
 
500
 
// array data layouts
501
 
// standard array layout  -----------------------------------------------------------------------------------------------
502
 
#define ALSTRING "Standard Array Layout"
 
699
#define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set) ]
 
700
#define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set) ]
 
701
 
 
702
// array handling  -----------------------------------------------------------------------------------------------
 
703
 
503
704
//#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
504
705
#define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
505
706
#define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
506
 
#define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set, l)*dTotalNum +(l)])
507
 
#define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set, l)*dTotalNum +(l)])
 
707
#define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set, l)*dTotalNum +(l)])
 
708
#define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set, l)*dTotalNum +(l)])
508
709
 
509
710
#define QCELLSTEP dTotalNum
510
711
#define _RACPNT(level, ii,ij,ik, is )  &QCELL(level,ii,ij,ik,is,0)
555
756
#define dNW 6
556
757
#define dSE 7
557
758
#define dSW 8
558
 
#define LBM_DFNUM 9
559
759
#else
560
760
// direction indices
561
761
#define dC 0
577
777
#define dEB 16
578
778
#define dWT 17
579
779
#define dWB 18
580
 
#define LBM_DFNUM 19
581
780
#endif
582
781
//? #define dWB 18
583
782
 
584
783
// default init for dFlux values
585
 
#define FLUX_INIT 0.5f * (float)(D::cDfNum)
 
784
#define FLUX_INIT 0.5f * (float)(this->cDfNum)
586
785
 
587
786
// only for non DF dir handling!
588
787
#define dNET 19
617
816
 
618
817
 
619
818
 
 
819
/******************************************************************************/
 
820
/*! equilibrium functions */
 
821
/******************************************************************************/
 
822
 
 
823
/*! calculate length of velocity vector */
 
824
inline LbmFloat LbmFsgrSolver::getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) {
 
825
        return ((ux)*dfDvecX[l]+(uy)*dfDvecY[l]+(uz)*dfDvecZ[l]);
 
826
};
 
827
 
 
828
/*! calculate equilibrium DF for given values */
 
829
inline LbmFloat LbmFsgrSolver::getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz) {
 
830
#if FSGR_STRICT_DEBUG==1
 
831
        if((l<0)||(l>LBM_DFNUM)) { errFatal("LbmFsgrSolver::getCollideEq","Invalid DFEQ call "<<l, SIMWORLD_PANIC ); /* no access to mPanic here */     }
 
832
#endif // FSGR_STRICT_DEBUG==1
 
833
        LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 
 
834
        return( dfLength[l] *( 
 
835
                                + rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
 
836
                                + 3.0 *tmp 
 
837
                                + 9.0/2.0 *(tmp*tmp) )
 
838
                        );
 
839
};
 
840
 
620
841
/*****************************************************************************/
621
842
/* init a given cell with flag, density, mass and equilibrium dist. funcs */
622
843
 
623
 
template<class D>
624
 
void LbmFsgrSolver<D>::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
 
844
void LbmFsgrSolver::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
625
845
        CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask;
626
846
        RFLAG(level,xx,yy,zz,set) = newflag | pers;
627
847
}
628
848
 
629
 
template<class D>
630
849
void 
631
 
LbmFsgrSolver<D>::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
 
850
LbmFsgrSolver::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
632
851
  /* init eq. dist funcs */
633
852
        LbmFloat *ecel;
634
853
        int workSet = mLevel[level].setCurr;
635
854
 
636
855
        ecel = RACPNT(level, i,j,k, workSet);
637
 
        FORDF0 { RAC(ecel, l) = D::dfEquil[l] * rho; }
 
856
        FORDF0 { RAC(ecel, l) = this->dfEquil[l] * rho; }
638
857
        RAC(ecel, dMass) = mass;
639
858
        RAC(ecel, dFfrac) = mass/rho;
640
859
        RAC(ecel, dFlux) = FLUX_INIT;
646
865
        return;
647
866
}
648
867
 
649
 
template<class D>
650
868
void 
651
 
LbmFsgrSolver<D>::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
 
869
LbmFsgrSolver::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
652
870
        LbmFloat *ecel;
653
871
        int workSet = mLevel[level].setCurr;
654
872
 
655
873
        ecel = RACPNT(level, i,j,k, workSet);
656
 
        FORDF0 { RAC(ecel, l) = D::getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
 
874
        FORDF0 { RAC(ecel, l) = getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
657
875
        RAC(ecel, dMass) = mass;
658
876
        RAC(ecel, dFfrac) = mass/rho;
659
877
        RAC(ecel, dFlux) = FLUX_INIT;
665
883
        return;
666
884
}
667
885
 
668
 
template<class D>
669
 
int LbmFsgrSolver<D>::getForZMinBnd() { 
 
886
int LbmFsgrSolver::getForZMinBnd() { 
670
887
        return 0; 
671
888
}
672
 
template<class D>
673
 
int LbmFsgrSolver<D>::getForZMin1()   { 
674
 
        if(D::cDimension==2) return 0;
 
889
int LbmFsgrSolver::getForZMin1()   { 
 
890
        if(LBMDIM==2) return 0;
675
891
        return 1; 
676
892
}
677
893
 
678
 
template<class D>
679
 
int LbmFsgrSolver<D>::getForZMaxBnd(int lev) { 
680
 
        if(D::cDimension==2) return 1;
 
894
int LbmFsgrSolver::getForZMaxBnd(int lev) { 
 
895
        if(LBMDIM==2) return 1;
681
896
        return mLevel[lev].lSizez -0;
682
897
}
683
 
template<class D>
684
 
int LbmFsgrSolver<D>::getForZMax1(int lev)   { 
685
 
        if(D::cDimension==2) return 1;
 
898
int LbmFsgrSolver::getForZMax1(int lev)   { 
 
899
        if(LBMDIM==2) return 1;
686
900
        return mLevel[lev].lSizez -1;
687
901
}
688
902