This version now includes the fluid control sources, however the Blender
[blender-staging.git] / intern / elbeem / intern / solver_control.h
1 /******************************************************************************
2  *
3  * El'Beem - the visual lattice boltzmann freesurface simulator
4  * All code distributed as part of El'Beem is covered by the version 2 of the 
5  * GNU General Public License. See the file COPYING for details.
6  * Copyright 2003-2006 Nils Thuerey
7  *
8  * testing extensions
9  *
10  *****************************************************************************/
11
12
13 #ifndef LBM_TESTCLASS_H
14 #define LBM_TESTCLASS_H
15
16 //class IsoSurface;
17 class ParticleObject;
18 class ControlParticles;
19 class ControlForces;
20
21 //#define NUMGRIDS 2
22 //#define MAXNUMSWS 10
23
24 // farfield modes
25 #define FARF_3DONLY -1
26 #define FARF_BOTH    0
27 #define FARF_SWEONLY 1
28 // dont reuse 3d vars/init
29 #define FARF_SEPSWE  2
30
31 // relaxation macros for solver_relax.h
32 #if LBM_INCLUDE_CONTROL!=1
33
34 // defined in relax.h
35
36 #else // LBM_INCLUDE_TESTSOLVERS!=1
37
38 // WARNING has to match controlparts.h
39 #define CPF_ENTRIES     12
40 #define CPF_FORCE       0
41 #define CPF_VELWEIGHT   3
42 #define CPF_VELOCITY    4
43 #define CPF_FORCEWEIGHT 7
44 #define CPF_MINCPDIST   8
45 #define CPF_MINCPDIR    9
46
47 #include "controlparticles.h"
48
49 // get force entry, set=0 is unused anyway
50 #define LBMGET_FORCE(lev, i,j,k)  mpControl->mCpForces[lev][ (LBMGI(lev,i,j,k,0)) ]
51
52 // debug mods off...
53 // same as in src/solver_relax.h!
54 #define __PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
55         ux += (grav)[0]; \
56         uy += (grav)[1]; \
57         uz += (grav)[2]; 
58
59 //void testMaxdmod(int i, int j,int k, LbmFloat &ux,LbmFloat &uy,LbmFloat &uz,ControlForces &ff);
60 #if LBMDIM==3
61 #define MAXDGRAV \
62                         if(myforce->forceMaxd[0]*ux+myforce->forceMaxd[1]*uy<LBM_EPSILON) { \
63                                 ux = v2w*myforce->forceVel[0]+ v2wi*ux;  \
64                                 uy = v2w*myforce->forceVel[1]+ v2wi*uy; }  \
65                         /* movement inverse to g? */ \
66                         if((uz>LBM_EPSILON)&&(uz>myforce->forceVel[2])) { \
67                                 uz = v2w*myforce->forceVel[2]+ v2wi*uz; } 
68 #else // LBMDIM==3
69 #define MAXDGRAV \
70                         if(myforce->forceMaxd[0]*ux<LBM_EPSILON) { \
71                                 ux = v2w*myforce->forceVel[0]+ v2wi*ux; } \
72                         /* movement inverse to g? */ \
73                         if((uy>LBM_EPSILON)&&(uy>myforce->forceVel[1])) { \
74                                 uy = v2w*myforce->forceVel[1]+ v2wi*uy; } 
75 #endif // LBMDIM==3
76
77 // debug modifications of collide vars (testing)
78 // requires: lev,i,j,k
79 #define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
80         LbmFloat attforce = 1.; \
81         if(this->mTForceStrength>0.) { \
82                 ControlForces* myforce = &LBMGET_FORCE(lev,i,j,k); \
83                 const LbmFloat vf = myforce->weightAtt;\
84                 const LbmFloat vw = myforce->weightVel;\
85                 if(vf!=0.) { attforce = MAX(0., 1.-vf);  /* TODO FIXME? use ABS(vf) for repulsion force? */ \
86                         ux += myforce->forceAtt[0]; \
87                         uy += myforce->forceAtt[1]; \
88                         uz += myforce->forceAtt[2]; \
89                         \
90                 } else if(( myforce->maxDistance>0.) && ( myforce->maxDistance<CPF_MAXDINIT)) {\
91                         const LbmFloat v2w = mpControl->mCons[0]->mCparts->getInfluenceMaxdist() * \
92                                 (myforce->maxDistance-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()) / (mpControl->mCons[0]->mCparts->getRadiusMaxd()-mpControl->mCons[0]->mCparts->getRadiusMinMaxd()); \
93                         const LbmFloat v2wi = 1.-v2w; \
94                         if(v2w>0.){ MAXDGRAV; \
95                                 /* errMsg("ERRMDTT","at "<<PRINT_IJK<<" maxd="<<myforce->maxDistance<<", newu"<<PRINT_VEC(ux,uy,uz)<<", org"<<PRINT_VEC(oux,ouy,ouz)<<", fv"<<myforce->forceVel<<" " );  */ \
96                         }\
97                 } \
98                 if(vw>0.) { \
99                         const LbmFloat vwi = 1.-vw;\
100                         const LbmFloat vwd = mpControl->mDiffVelCon;\
101                         ux += vw*(myforce->forceVel[0]-myforce->compAv[0] + vwd*(myforce->compAv[0]-ux) ); \
102                         uy += vw*(myforce->forceVel[1]-myforce->compAv[1] + vwd*(myforce->compAv[1]-uy) ); \
103                         uz += vw*(myforce->forceVel[2]-myforce->compAv[2] + vwd*(myforce->compAv[2]-uz) ); \
104                         /*  TODO test!? modify smooth vel by influence of force for each lbm step, to account for force update only each N steps */ \
105                         myforce->compAv = (myforce->forceVel*vw+ myforce->compAv*vwi); \
106                 } \
107         } \
108         ux += (grav)[0]*attforce; \
109         uy += (grav)[1]*attforce; \
110         uz += (grav)[2]*attforce; \
111         /* end PRECOLLIDE_MODS */
112
113 #define TEST_IF_CHECK \
114                 if((!iffilled)&&(LBMGET_FORCE(lev,i,j,k).weightAtt!=0.)) { \
115                         errMsg("TESTIFFILL"," at "<<PRINT_IJK<<" "<<mass<<" "<<rho); \
116                         iffilled = true; \
117                         if(mass<rho*1.0) mass = rho*1.0; myfrac = 1.0; \
118                 }
119
120 #endif // LBM_INCLUDE_TESTSOLVERS!=1
121
122
123 // a single set of control particles and params
124 class LbmControlSet {
125         public:
126                 LbmControlSet();
127                 ~LbmControlSet();
128                 void initCparts();
129
130                 // control particles
131                 ControlParticles *mCparts; 
132                 // control particle overall motion (for easier manual generation)
133                 ControlParticles *mCpmotion;
134                 // cp data file
135                 string mContrPartFile;
136                 string mCpmotionFile;
137                 // cp debug displau
138                 LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
139
140                 // params
141                 AnimChannel<float> mcForceAtt;
142                 AnimChannel<float> mcForceVel;
143                 AnimChannel<float> mcForceMaxd;
144
145                 AnimChannel<float> mcRadiusAtt;
146                 AnimChannel<float> mcRadiusVel;
147                 AnimChannel<float> mcRadiusMind;
148                 AnimChannel<float> mcRadiusMaxd;
149
150                 AnimChannel<ntlVec3f> mcCpScale;
151                 AnimChannel<ntlVec3f> mcCpOffset;
152 };
153                 
154
155
156 // main control data storage
157 class LbmControlData 
158 {
159         public:
160                 LbmControlData();
161                 virtual ~LbmControlData();
162
163                 // control data
164
165                 // contorl params
166                 void parseControldataAttrList(AttributeList *attr);
167
168                 // control strength, set for solver interface
169                 LbmFloat mSetForceStrength;
170                 // cp vars
171                 std::vector<LbmControlSet*> mCons;
172                 // update interval
173                 int mCpUpdateInterval;
174                 // output
175                 string mCpOutfile;
176                 // control particle precomputed influence
177                 std::vector< std::vector<ControlForces> > mCpForces;
178                 std::vector<ControlForces> mCpKernel;
179                 std::vector<ControlForces> mMdKernel;
180                 // activate differential velcon
181                 LbmFloat mDiffVelCon;
182
183                 // cp debug displau
184                 LbmFloat mDebugCpscale, mDebugVelScale, mDebugCompavScale, mDebugAttScale, mDebugMaxdScale, mDebugAvgVelScale;
185 };
186
187 #endif // LBM_TESTCLASS_H