merge from trunk #37722
[blender-staging.git] / intern / elbeem / intern / solver_interface.h
1 /** \file elbeem/intern/solver_interface.h
2  *  \ingroup elbeem
3  */
4 /******************************************************************************
5  *
6  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
7  * All code distributed as part of El'Beem is covered by the version 2 of the 
8  * GNU General Public License. See the file COPYING for details.
9  * Copyright 2003-2006 Nils Thuerey
10  *
11  * Header for Combined 2D/3D Lattice Boltzmann Interface Class
12  * 
13  *****************************************************************************/
14 #ifndef LBMINTERFACE_H
15 #define LBMINTERFACE_H
16
17 //! include gui support?
18 #ifndef NOGUI
19 #define LBM_USE_GUI      1
20 #else
21 #define LBM_USE_GUI      0
22 #endif
23
24 #if LBM_USE_GUI==1
25 #define USE_GLUTILITIES
26 // for debug display
27 //#include <GL/gl.h>
28 #include "../gui/guifuncs.h"
29 #endif
30
31 #include <sstream>
32 #include "utilities.h"
33 #include "ntl_bsptree.h"
34 #include "ntl_geometryobject.h"
35 #include "parametrizer.h"
36 #include "attributes.h"
37 #include "isosurface.h"
38 class ParticleTracer;
39 class ParticleObject;
40
41 // use which fp-precision for LBM? 1=float, 2=double
42 #ifdef PRECISION_LBM_SINGLE
43 #define LBM_PRECISION 1
44 #else
45 #ifdef PRECISION_LBM_DOUBLE
46 #define LBM_PRECISION 2
47 #else
48 // default to floats
49 #define LBM_PRECISION 1
50 #endif
51 #endif
52
53 #if LBM_PRECISION==1
54 /* low precision for LBM solver */
55 typedef float LbmFloat;
56 typedef ntlVec3f LbmVec;
57 #define LBM_EPSILON (1e-5)
58 #else
59 /* standard precision for LBM solver */
60 typedef double LbmFloat;
61 typedef ntlVec3d LbmVec;
62 #define LBM_EPSILON (1e-10)
63 #endif
64
65 // long integer, needed e.g. for memory calculations
66 #ifndef USE_MSVC6FIXES
67 #define LONGINT long long int
68 #else
69 #define LONGINT _int64
70 #endif
71
72
73 // default to 3dim
74 #ifndef LBMDIM
75 #define LBMDIM 3
76 #endif // LBMDIM
77
78 #if LBMDIM==2
79 #define LBM_DFNUM 9
80 #else
81 #define LBM_DFNUM 19
82 #endif
83
84 // conversions (lbm and parametrizer)
85 template<class T> inline LbmVec     vec2L(T v) { return LbmVec(v[0],v[1],v[2]); }
86 template<class T> inline ParamVec   vec2P(T v) { return ParamVec(v[0],v[1],v[2]); }
87
88 template<class Scalar> class ntlMatrix4x4;
89
90
91 // bubble id type
92 typedef int BubbleId;
93
94 // basic cell type distinctions
95 #define CFUnused              (1<< 0)
96 #define CFEmpty               (1<< 1)
97 #define CFBnd                 (1<< 2)
98 #define CFMbndInflow          (1<< 3)
99 #define CFMbndOutflow         (1<< 4)
100 #define CFFluid               (1<< 5)
101 #define CFInter               (1<< 6)
102 // additional for fluid (needed separately for adaptive grids)
103 #define CFNoBndFluid          (1<< 7)
104 #define CFNoDelete            (1<< 8)
105
106 // additional bnd add flags
107 #define CFBndNoslip           (1<< 9)
108 #define CFBndFreeslip         (1<<10)
109 #define CFBndPartslip         (1<<11)
110 #define CFBndMoving           (1<<12)
111
112 // additional for fluid/interface
113 // force symmetry for flag reinit 
114 #define CFNoInterpolSrc       (1<<13) 
115 #define CFNoNbFluid           (1<<14)
116 #define CFNoNbEmpty           (1<<15)
117         
118 // cell treated normally on coarser grids
119 #define CFGrNorm              (1<<16)
120 #define CFGrCoarseInited      (1<<17)
121
122 // (the following values shouldnt overlap to ensure
123 // proper coarsening)
124 // border cells to be interpolated from finer grid
125 #define CFGrFromFine          (1<<18)
126 // 32k aux border marker 
127 #define CFGrToFine            (1<<19)
128 // also needed on finest level
129 #define CFGrFromCoarse        (1<<20)
130 // additional refinement tags (coarse grids only?)
131 // */
132
133 // above 24 is used to encode in/outflow object type
134 #define CFPersistMask (0xFF000000 | CFMbndInflow | CFMbndOutflow)
135 #define CFNoPersistMask (~CFPersistMask)
136
137
138 // nk
139 #define CFInvalid             (CellFlagType)(1<<31)
140
141 // use 32bit flag types
142 //#ifdef __x86_64__
143  //typedef int cfINT32;
144 //#else
145  //typedef long cfINT32;
146 //#endif // defined (_IA64)  
147 //#define CellFlagType cfINT32
148 #define CellFlagType int
149 #define CellFlagTypeSize 4
150
151
152 // aux. field indices (same for 2d)
153 #define dFfrac 19
154 #define dMass 20
155 #define dFlux 21
156 // max. no. of cell values for 3d
157 #define dTotalNum 22
158
159
160 /*****************************************************************************/
161 /*! a single lbm cell */
162 /*  the template is only needed for 
163  *  dimension dependend constants e.g. 
164  *  number of df's in model */
165 class LbmCellContents {
166         public:
167                 LbmFloat     df[ 27 ]; // be on the safe side here...
168         LbmFloat     rho;
169                 LbmVec       vel;
170         LbmFloat     mass;
171                 CellFlagType flag;
172                 BubbleId     bubble;
173         LbmFloat     ffrac;
174 };
175
176 /* struct for the coordinates of a cell in the grid */
177 typedef struct {
178   int x,y,z;
179         int flag; // special handling?
180 } LbmPoint;
181
182 /* struct for the coordinates of a cell in the grid */
183 typedef struct {
184         char active;            // bubble in use, oder may be overwritten?
185   LbmFloat volume;          // volume of this bubble (0 vor atmosphere)
186         LbmFloat mass;            // "mass" of bubble 
187         int i,j,k;              // index of a cell in the bubble
188 } LbmBubble;
189
190
191
192
193 //! choose which data to display
194 #define FLUIDDISPINVALID    0
195 #define FLUIDDISPNothing    1
196 #define FLUIDDISPCelltypes  2
197 #define FLUIDDISPVelocities 3
198 #define FLUIDDISPCellfills  4
199 #define FLUIDDISPDensity    5
200 #define FLUIDDISPGrid       6
201 #define FLUIDDISPSurface    7
202
203
204
205 /*****************************************************************************/
206 //! cell identifier interface
207 class CellIdentifierInterface {
208         public:
209                 //! reset constructor
210                 CellIdentifierInterface():mEnd(false) { };
211                 //! virtual destructor
212                 virtual ~CellIdentifierInterface() {};
213
214                 //! return node as string (with some basic info)
215                 virtual string getAsString() = 0;
216
217                 //! compare cids
218                 virtual bool equal(CellIdentifierInterface* other) = 0;
219
220                 //! set/get end flag for grid traversal (not needed for marked cells)
221                 inline void setEnd(bool set){ mEnd = set;  }
222                 inline bool getEnd( )       { return mEnd; }
223
224                 //! has the grid been traversed?
225                 bool mEnd;
226
227 };
228
229
230
231 /*****************************************************************************/
232 /*! class defining abstract function interface */
233 /*  has to provide iterating functionality */
234 class LbmSolverInterface 
235 {
236         public:
237                 //! Constructor 
238                 LbmSolverInterface();
239                 //! Destructor 
240                 virtual ~LbmSolverInterface();
241                 //! id string of solver
242                 virtual string getIdString() = 0;
243
244                 // multi step solver init
245                 /*! finish the init with config file values (allocate arrays...) */
246                 virtual bool initializeSolverMemory() =0;
247                 /*! init solver arrays */
248                 virtual bool initializeSolverGrids() =0;
249                 /*! prepare actual simulation start, setup viz etc */
250                 virtual bool initializeSolverPostinit() =0;
251                 
252                 /*! notify object that dump is in progress (e.g. for field dump) */
253                 virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) = 0;
254
255                 /*! parse a boundary flag string */
256                 CellFlagType readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed);
257                 /*! parse standard attributes */
258                 void parseStdAttrList();
259                 /*! initilize variables fom attribute list (should at least call std parse) */
260                 virtual void parseAttrList() = 0;
261
262                 virtual void step() = 0;
263                 virtual void prepareVisualization() { /* by default off */ };
264
265                 /*! particle handling */
266                 virtual int initParticles() = 0;
267                 virtual void advanceParticles() = 0;
268                 /*! get surface object (NULL if no surface) */
269                 IsoSurface* getSurfaceGeoObj() { return mpIso; }
270
271                 /*! debug object display */
272                 virtual vector<ntlGeometryObject*> getDebugObjects() { vector<ntlGeometryObject*> empty(0); return empty; }
273
274                 /* surface generation settings */
275                 virtual void setSurfGenSettings(short value) = 0;
276
277 #if LBM_USE_GUI==1
278                 /*! show simulation info */
279                 virtual void debugDisplay(int) = 0;
280 #endif
281
282                 /*! init tree for certain geometry init */
283                 void initGeoTree();
284                 /*! destroy tree etc. when geometry init done */
285                 void freeGeoTree();
286                 /*! check for a certain flag type at position org (needed for e.g. quadtree refinement) */
287                 bool geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance,int shootDir=0);
288                 bool geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, 
289                                 const gfxReal halfCellsize, bool &thinHit, bool recurse);
290                 /*! set render globals, for scene/tree access */
291                 void setRenderGlobals(ntlRenderGlobals *glob) { mpGlob = glob; };
292                 /*! get max. velocity of all objects to initialize as fluid regions, and of all moving objects */
293                 ntlVec3Gfx getGeoMaxMovementVelocity(LbmFloat simtime, LbmFloat stepsize);
294
295                 /* rt interface functions */
296                 unsigned int getIsoVertexCount()  { return mpIso->getIsoVertexCount(); }
297                 unsigned int getIsoIndexCount()   { return mpIso->getIsoIndexCount(); }
298                 char* getIsoVertexArray()         { return mpIso->getIsoVertexArray(); }
299                 unsigned int *getIsoIndexArray()  { return mpIso->getIsoIndexArray(); }
300                 void triangulateSurface()         { mpIso->triangulate(); }
301
302                 /* access functions */
303
304                 /*! return grid sizes */
305                 int getSizeX( void ) { return mSizex; }
306                 int getSizeY( void ) { return mSizey; }
307                 int getSizeZ( void ) { return mSizez; }
308                 /*! return grid sizes */
309                 void setSizeX( int ns ) { mSizex = ns; }
310                 void setSizeY( int ns ) { mSizey = ns; }
311                 void setSizeZ( int ns ) { mSizez = ns; }
312                 /*! access fluid only simulation flag */
313                 void setAllfluid(bool set) { mAllfluid=set; }
314                 bool getAllfluid()         { return mAllfluid; }
315
316                 /*! set attr list pointer */
317                 void setAttrList(AttributeList *set) { mpSifAttrs = set; }
318                 /*! Returns the attribute list pointer */
319                 inline AttributeList *getAttributeList() { return mpSifAttrs; }
320                 /*! set sws attr list pointer */
321                 void setSwsAttrList(AttributeList *set) { mpSifSwsAttrs = set; }
322                 inline AttributeList *getSwsAttributeList() { return mpSifSwsAttrs; }
323
324                 /*! set parametrizer pointer */
325                 inline void setParametrizer(Parametrizer *set) { mpParam = set; }
326                 /*! get parametrizer pointer */
327                 inline Parametrizer *getParametrizer() { return mpParam; }
328                 /*! get/set particle pointer */
329                 inline void setParticleTracer(ParticleTracer *set) { mpParticles = set; }
330                 inline ParticleTracer *getParticleTracer() { return mpParticles; }
331
332                 /*! set density gradient init from e.g. init test cases */
333                 inline void setInitDensityGradient(bool set) { mInitDensityGradient = set; }
334
335                 /*! access geometry start vector */
336                 inline void setGeoStart(ntlVec3Gfx set) { mvGeoStart = set; }
337                 inline ntlVec3Gfx getGeoStart() const   { return mvGeoStart; }
338
339                 /*! access geometry end vector */
340                 inline void setGeoEnd(ntlVec3Gfx set)   { mvGeoEnd = set; }
341                 inline ntlVec3Gfx getGeoEnd() const     { return mvGeoEnd; }
342
343                 /*! access geo init vars */
344                 inline void setLbmInitId(int set)       { mLbmInitId = set; }
345                 inline int getLbmInitId() const { return mLbmInitId; }
346
347                 /*! init domain transformation matrix from float array */
348                 void initDomainTrafo(float *mat);
349                 /*! get domain transformation matrix to have object centered fluid vertices */
350                 inline ntlMatrix4x4<gfxReal> *getDomainTrafo() { return mpSimTrafo; }
351
352                 /*! access name string */
353                 inline void setName(string set) { mName = set; }
354                 inline string getName() const   { return mName; }
355
356                 /*! access string for node info debugging output */
357                 inline string getNodeInfoString() const { return mNodeInfoString; }
358
359                 /*! get panic flag */
360                 inline bool getPanic() { return mPanic; }
361
362                 //! set silent mode?
363                 inline void setSilent(bool set){ mSilent = set; }
364
365                 //! set amount of surface/normal smoothing
366                 inline void setSmoothing(float setss,float setns){ mSmoothSurface=setss; mSmoothNormals=setns; }
367                 //! set amount of iso subdivisions
368                 inline void setIsoSubdivs(int s){ mIsoSubdivs=s; }
369                 //! set desired refinement
370                 inline void setPreviewSize(int set){ mOutputSurfacePreview = set; }
371                 //! set desired refinement
372                 inline void setRefinementDesired(int set){ mRefinementDesired = set; }
373
374                 //! set/get dump velocities flag
375                 inline void setDumpVelocities(bool set) { mDumpVelocities = set; }
376                 inline bool getDumpVelocities() const   { return mDumpVelocities; }
377
378                 //! set/get particle generation prob.
379                 inline void setGenerateParticles(LbmFloat set)  { mPartGenProb = set; }
380                 inline LbmFloat getGenerateParticles() const    { return mPartGenProb; }
381
382                 //! set/get dump velocities flag
383                 inline void setDomainBound(string set)  { mDomainBound = set; }
384                 inline string getDomainBound() const    { return mDomainBound; }
385                 //! set/get dump velocities flag
386                 inline void setDomainPartSlip(LbmFloat set)     { mDomainPartSlipValue = set; }
387                 inline LbmFloat getDomainPartSlip() const       { return mDomainPartSlipValue; }
388                 //! set/get far field size
389                 inline void setFarFieldSize(LbmFloat set)       { mFarFieldSize = set; }
390                 inline LbmFloat getFarFieldSize() const { return mFarFieldSize; }
391                 //! set/get cp stage
392                 inline void setCpStage(int set) { mCppfStage = set; }
393                 inline int getCpStage() const   { return mCppfStage; }
394                 //! set/get dump modes
395                 inline void setDumpRawText(bool set)    { mDumpRawText = set; }
396                 inline bool getDumpRawText() const      { return mDumpRawText; }
397                 inline void setDumpRawBinary(bool set)  { mDumpRawBinary = set; }
398                 inline bool getDumpRawBinary() const    { return mDumpRawBinary; }
399                 inline void setDumpRawBinaryZip(bool set)       { mDumpRawBinaryZip = set; }
400                 inline bool getDumpRawBinaryZip() const { return mDumpRawBinaryZip; }
401                 //! set/get debug vel scale
402                 inline void setDebugVelScale(LbmFloat set)      { mDebugVelScale = set; }
403                 inline LbmFloat getDebugVelScale() const        { return mDebugVelScale; }
404
405                 // cell iterator interface
406                 
407                 // cell id type
408                 typedef CellIdentifierInterface* CellIdentifier;
409
410                 //! cell iteration methods
411                 virtual CellIdentifierInterface* getFirstCell( ) = 0;
412                 virtual void advanceCell( CellIdentifierInterface* ) = 0;
413                 virtual bool noEndCell( CellIdentifierInterface* ) = 0;
414                 //! clean up iteration, this should be called, when the iteration is not completely finished
415                 virtual void deleteCellIterator( CellIdentifierInterface** ) = 0;
416
417                 //! find cell at a given position (returns NULL if not in domain)
418                 virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos ) = 0;
419
420                 //! return node information
421                 virtual int        getCellSet      ( CellIdentifierInterface* ) = 0;
422                 virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* ) = 0;
423                 virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* ) = 0;
424                 virtual int        getCellLevel    ( CellIdentifierInterface* ) = 0;
425                 virtual LbmFloat   getCellDensity  ( CellIdentifierInterface*,int ) = 0;
426                 virtual LbmVec     getCellVelocity ( CellIdentifierInterface*,int ) = 0;
427                 /*! get equilibrium distribution functions */
428                 virtual LbmFloat   getEquilDf      ( int ) = 0;
429                 /*! redundant cell functions */
430                 virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir) = 0;
431                 virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set) = 0;
432                 virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set) = 0;
433                 virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set) = 0;
434
435                 /*! get velocity directly from position */
436                 virtual ntlVec3Gfx getVelocityAt(float x, float y, float z) = 0;
437
438                 // gui/output debugging functions
439 #if LBM_USE_GUI==1
440                 virtual void debugDisplayNode(int dispset, CellIdentifier cell ) = 0;
441                 virtual void lbmDebugDisplay(int dispset) = 0;
442                 virtual void lbmMarkedCellDisplay() = 0;
443 #endif // LBM_USE_GUI==1
444                 virtual void debugPrintNodeInfo(CellIdentifier cell, int forceSet=-1) = 0;
445
446                 // debugging cell marker functions
447
448                 //! add cell to mMarkedCells list
449                 void addCellToMarkedList( CellIdentifierInterface *cid );
450                 //! marked cell iteration methods
451                 CellIdentifierInterface* markedGetFirstCell( );
452                 CellIdentifierInterface* markedAdvanceCell();
453                 void markedClearList();
454
455
456         protected:
457
458                 /*! abort simulation on error... */
459                 bool mPanic;
460
461
462                 /*! Size of the array in x,y,z direction */
463                 int mSizex, mSizey, mSizez;
464                 /*! only fluid in sim? */
465                 bool mAllfluid;
466
467
468                 /*! step counter */
469                 int mStepCnt;
470
471                 /*! mass change from one step to the next, for extreme cases fix globally */
472                 LbmFloat mFixMass;
473
474                 // deprecated param vars
475                 /*! omega for lbm */
476                 LbmFloat mOmega;
477                 /*! gravity strength in neg. z direction */
478                 LbmVec mGravity;
479                 /*! Surface tension of the fluid */
480                 LbmFloat mSurfaceTension;
481
482
483                 /* boundary inits */
484                 CellFlagType mBoundaryEast, mBoundaryWest, 
485                   mBoundaryNorth, mBoundarySouth, 
486                   mBoundaryTop, mBoundaryBottom;
487
488                 /*! initialization from config file done? */
489                 int mInitDone;
490
491                 /*! init density gradient? */
492                 bool mInitDensityGradient;
493
494                 /*! pointer to the attribute list, only pointer to obj attrs */
495                 AttributeList *mpSifAttrs;
496                 AttributeList *mpSifSwsAttrs;
497
498                 /*! get parameters from this parametrize in finishInit */
499                 Parametrizer *mpParam;
500                 //! store particle tracer
501                 ParticleTracer *mpParticles;
502
503                 /*! number of particles lost so far */
504                 int mNumParticlesLost;
505                 /*! number of particles lost so far */
506                 int mNumInvalidDfs;
507                 /*! no of filled/emptied cells per time step */
508                 int mNumFilledCells, mNumEmptiedCells;
509                 /*! counter number of used cells for performance */
510                 int mNumUsedCells;
511                 /*! MLSUPS counter */
512                 LbmFloat mMLSUPS;
513                 /*! debug - velocity output scaling factor */
514                 LbmFloat mDebugVelScale;
515                 /*! string for node info debugging output */
516                 string mNodeInfoString;
517
518                 // geo init vars
519                 // TODO deprecate SimulationObject vars
520
521                 /*! for display - start and end vectors for geometry */
522                 ntlVec3Gfx mvGeoStart, mvGeoEnd;
523                 //! domain vertex trafos
524                 ntlMatrix4x4<gfxReal> *mpSimTrafo;
525
526                 /*! perform accurate geometry init? */
527                 bool mAccurateGeoinit;
528
529                 /*! name of this lbm object (for debug output) */
530                 string mName;
531
532                 //! Mcubes object for surface reconstruction 
533                 IsoSurface *mpIso;
534                 /*! isolevel value for marching cubes surface reconstruction */
535                 LbmFloat mIsoValue;
536
537                 //! debug output?
538                 bool mSilent;
539
540                 /*! geometry init id, passed from ntl_geomclass */
541                 int mLbmInitId;
542                 /*! tree object for geomerty initialization */
543                 ntlTree *mpGiTree;
544                 /*! object vector for geo init */
545                 vector<ntlGeometryObject*> *mpGiObjects;
546                 /*! inside which objects? */
547                 vector<int> mGiObjInside;
548                 /*! inside which objects? */
549                 vector<gfxReal> mGiObjDistance;
550                 vector<gfxReal> mGiObjSecondDist;
551                 /*! remember globals */
552                 ntlRenderGlobals *mpGlob;
553                 
554                 //! use refinement/coarsening?
555                 int mRefinementDesired;
556
557                 //! output surface preview? if >0 yes, and use as reduzed size 
558                 int mOutputSurfacePreview;
559                 LbmFloat mPreviewFactor;
560
561                 /*! enable surface and normals smoothing? */
562                 float mSmoothSurface;
563                 float mSmoothNormals;
564                 /*! isosurface subdivisions */
565                 int mIsoSubdivs;
566
567                 //! particle generation probability
568                 LbmFloat mPartGenProb;
569
570                 //! dump velocities?
571                 bool mDumpVelocities;
572
573                 // list for marked cells
574                 vector<CellIdentifierInterface *> mMarkedCells;
575                 int mMarkedCellIndex;
576
577                 //! domain boundary free/no slip type
578                 string mDomainBound;
579                 //! part slip value for domain
580                 LbmFloat mDomainPartSlipValue;
581
582                 // size of far field area
583                 LbmFloat mFarFieldSize;
584                 // amount of drop mass to subtract
585                 LbmFloat mPartDropMassSub;
586                 // use physical drop model for particles?
587                 bool mPartUsePhysModel;
588
589                 //! test vars
590                 // strength of applied force
591                 LbmFloat mTForceStrength;
592                 int mCppfStage;
593
594                 //! dumping modes
595                 bool mDumpRawText;
596                 bool mDumpRawBinary;
597                 bool mDumpRawBinaryZip;
598 };
599
600
601 // helper function to create consistent grid resolutions
602 void initGridSizes(int &mSizex, int &mSizey, int &mSizez,
603                 ntlVec3Gfx &mvGeoStart, ntlVec3Gfx &mvGeoEnd, 
604                 int mMaxRefine, bool parallel);
605 // return the amount of memory required in total (reqret)
606 // and for the finest grid only (reqretFine, can be NULL)
607 void calculateMemreqEstimate(int resx,int resy,int resz, int refine,
608                 float farfieldsize, double *reqret, double *reqretFine, string *reqstr);
609
610 //! helper function to convert flag to string (for debuggin)
611 string convertCellFlagType2String( CellFlagType flag );
612 string convertSingleFlag2String(CellFlagType cflag);
613
614 #endif // LBMINTERFACE_H