doxygen: intern/elbeem tagged
[blender.git] / intern / elbeem / intern / solver_control.cpp
1 /** \file elbeem/intern/solver_control.cpp
2  *  \ingroup elbeem
3  */
4 /******************************************************************************
5  *
6  * El'Beem - the visual lattice boltzmann freesurface simulator
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  *
10  * Copyright 2003-2008 Nils Thuerey 
11  *
12  * control extensions
13  *
14  *****************************************************************************/
15 #include "solver_class.h"
16 #include "solver_relax.h"
17 #include "particletracer.h"
18
19 #include "solver_control.h"
20
21 #include "controlparticles.h"
22
23 #include "elbeem.h"
24
25 #include "ntl_geometrymodel.h"
26
27 /******************************************************************************
28  * LbmControlData control set
29  *****************************************************************************/
30
31 LbmControlSet::LbmControlSet() :
32         mCparts(NULL), mCpmotion(NULL), mContrPartFile(""), mCpmotionFile(""),
33         mcForceAtt(0.), mcForceVel(0.), mcForceMaxd(0.), 
34         mcRadiusAtt(0.), mcRadiusVel(0.), mcRadiusMind(0.), mcRadiusMaxd(0.), 
35         mcCpScale(1.), mcCpOffset(0.)
36 {
37 }
38 LbmControlSet::~LbmControlSet() {
39         if(mCparts) delete mCparts;
40         if(mCpmotion) delete mCpmotion;
41 }
42 void LbmControlSet::initCparts() {
43         mCparts = new ControlParticles();
44         mCpmotion = new ControlParticles();
45 }
46
47
48
49 /******************************************************************************
50  * LbmControlData control
51  *****************************************************************************/
52
53 LbmControlData::LbmControlData() : 
54         mSetForceStrength(0.),
55         mCons(), 
56         mCpUpdateInterval(8), // DG: was 16 --> causes problems (big sphere after some time), unstable
57         mCpOutfile(""), 
58         mCpForces(), mCpKernel(), mMdKernel(),
59         mDiffVelCon(1.),
60         mDebugCpscale(0.), 
61         mDebugVelScale(0.), 
62         mDebugCompavScale(0.), 
63         mDebugAttScale(0.), 
64         mDebugMaxdScale(0.),
65         mDebugAvgVelScale(0.)
66 {
67 }
68
69 LbmControlData::~LbmControlData() 
70 {
71         while (!mCons.empty()) {
72                 delete mCons.back();  mCons.pop_back();
73         }
74 }
75
76
77 void LbmControlData::parseControldataAttrList(AttributeList *attr) {
78         // controlpart vars
79         mSetForceStrength = attr->readFloat("tforcestrength", mSetForceStrength,"LbmControlData", "mSetForceStrength", false);
80         //errMsg("tforcestrength set to "," "<<mSetForceStrength);
81         mCpUpdateInterval = attr->readInt("controlparticle_updateinterval", mCpUpdateInterval,"LbmControlData","mCpUpdateInterval", false);
82         // tracer output file
83         mCpOutfile = attr->readString("controlparticle_outfile",mCpOutfile,"LbmControlData","mCpOutfile", false);
84         if(getenv("ELBEEM_CPOUTFILE")) {
85                 string outfile(getenv("ELBEEM_CPOUTFILE"));
86                 mCpOutfile = outfile;
87                 debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPOUTFILE to set mCpOutfile to "<<outfile<<","<<mCpOutfile,7);
88         }
89
90         for(int cpii=0; cpii<10; cpii++) {
91                 string suffix("");
92                 //if(cpii>0)
93                 {  suffix = string("0"); suffix[0]+=cpii; }
94                 LbmControlSet *cset;
95                 cset = new LbmControlSet();
96                 cset->initCparts();
97
98                 cset->mContrPartFile = attr->readString("controlparticle"+suffix+"_file",cset->mContrPartFile,"LbmControlData","cset->mContrPartFile", false);
99                 if((cpii==0) && (getenv("ELBEEM_CPINFILE")) ) {
100                         string infile(getenv("ELBEEM_CPINFILE"));
101                         cset->mContrPartFile = infile;
102                         debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPINFILE to set mContrPartFile to "<<infile<<","<<cset->mContrPartFile,7);
103                 }
104
105                 LbmFloat cpvort=0.;
106                 cset->mcRadiusAtt =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusatt", 0., "LbmControlData","mcRadiusAtt" );
107                 cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
108                 cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
109                 cset->mCparts->setRadiusAtt(cset->mcRadiusAtt.get(0.));
110                 cset->mCparts->setRadiusVel(cset->mcRadiusVel.get(0.));
111
112                 // WARNING currently only for first set
113                 //if(cpii==0) {
114                 cset->mcForceAtt  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_attraction", 0. , "LbmControlData","cset->mcForceAtt", false);
115                 cset->mcForceVel  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_velocity",   0. , "LbmControlData","mcForceVel", false);
116                 cset->mcForceMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_maxdist",    0. , "LbmControlData","mcForceMaxd", false);
117                 cset->mCparts->setInfluenceAttraction(cset->mcForceAtt.get(0.) );
118                 // warning - stores temprorarily, value converted to dt dep. factor
119                 cset->mCparts->setInfluenceVelocity(cset->mcForceVel.get(0.) , 0.01 ); // dummy dt
120                 cset->mCparts->setInfluenceMaxdist(cset->mcForceMaxd.get(0.) );
121                 cpvort =  attr->readFloat("controlparticle"+suffix+"_vorticity",   cpvort, "LbmControlData","cpvort", false);
122                 cset->mCparts->setInfluenceTangential(cpvort);
123                         
124                 cset->mcRadiusMind =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmin", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMind", false);
125                 cset->mcRadiusMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmax", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMaxd", false);
126                 cset->mCparts->setRadiusMinMaxd(cset->mcRadiusMind.get(0.));
127                 cset->mCparts->setRadiusMaxd(cset->mcRadiusMaxd.get(0.));
128                 //}
129
130                 // now local...
131                 //LbmVec cpOffset(0.), cpScale(1.);
132                 LbmFloat cpTimescale = 1.;
133                 string cpMirroring("");
134
135                 //cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
136                 //cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
137                 cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
138                 cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
139                 cpTimescale =  attr->readFloat("controlparticle"+suffix+"_timescale",  cpTimescale, "LbmControlData","cpTimescale", false);
140                 cpMirroring =  attr->readString("controlparticle"+suffix+"_mirror",  cpMirroring, "LbmControlData","cpMirroring", false);
141
142                 LbmFloat cpsWidth = cset->mCparts->getCPSWith();
143                 cpsWidth =  attr->readFloat("controlparticle"+suffix+"_cpswidth",  cpsWidth, "LbmControlData","cpsWidth", false);
144                 LbmFloat cpsDt = cset->mCparts->getCPSTimestep();
145                 cpsDt =  attr->readFloat("controlparticle"+suffix+"_cpstimestep",  cpsDt, "LbmControlData","cpsDt", false);
146                 LbmFloat cpsTstart = cset->mCparts->getCPSTimeStart();
147                 cpsTstart =  attr->readFloat("controlparticle"+suffix+"_cpststart",  cpsTstart, "LbmControlData","cpsTstart", false);
148                 LbmFloat cpsTend = cset->mCparts->getCPSTimeEnd();
149                 cpsTend =  attr->readFloat("controlparticle"+suffix+"_cpstend",  cpsTend, "LbmControlData","cpsTend", false);
150                 LbmFloat cpsMvmfac = cset->mCparts->getCPSMvmWeightFac();
151                 cpsMvmfac =  attr->readFloat("controlparticle"+suffix+"_cpsmvmfac",  cpsMvmfac, "LbmControlData","cpsMvmfac", false);
152                 cset->mCparts->setCPSWith(cpsWidth);
153                 cset->mCparts->setCPSTimestep(cpsDt);
154                 cset->mCparts->setCPSTimeStart(cpsTstart);
155                 cset->mCparts->setCPSTimeEnd(cpsTend);
156                 cset->mCparts->setCPSMvmWeightFac(cpsMvmfac);
157
158                 cset->mCparts->setOffset( vec2L(cset->mcCpOffset.get(0.)) );
159                 cset->mCparts->setScale( vec2L(cset->mcCpScale.get(0.)) );
160                 cset->mCparts->setInitTimeScale( cpTimescale );
161                 cset->mCparts->setInitMirror( cpMirroring );
162
163                 int mDebugInit = 0;
164                 mDebugInit = attr->readInt("controlparticle"+suffix+"_debuginit", mDebugInit,"LbmControlData","mDebugInit", false);
165                 cset->mCparts->setDebugInit(mDebugInit);
166
167                 // motion particle settings
168                 LbmVec mcpOffset(0.), mcpScale(1.);
169                 LbmFloat mcpTimescale = 1.;
170                 string mcpMirroring("");
171
172                 cset->mCpmotionFile = attr->readString("cpmotion"+suffix+"_file",cset->mCpmotionFile,"LbmControlData","mCpmotionFile", false);
173                 mcpTimescale =  attr->readFloat("cpmotion"+suffix+"_timescale",  mcpTimescale, "LbmControlData","mcpTimescale", false);
174                 mcpMirroring =  attr->readString("cpmotion"+suffix+"_mirror",  mcpMirroring, "LbmControlData","mcpMirroring", false);
175                 mcpOffset = vec2L( attr->readVec3d("cpmotion"+suffix+"_offset", vec2P(mcpOffset),"LbmControlData","cpOffset", false) );
176                 mcpScale =  vec2L( attr->readVec3d("cpmotion"+suffix+"_scale",  vec2P(mcpScale), "LbmControlData","cpScale", false) );
177
178                 cset->mCpmotion->setOffset( vec2L(mcpOffset) );
179                 cset->mCpmotion->setScale( vec2L(mcpScale) );
180                 cset->mCpmotion->setInitTimeScale( mcpTimescale );
181                 cset->mCpmotion->setInitMirror( mcpMirroring );
182
183                 if(cset->mContrPartFile.length()>1) {
184                         errMsg("LbmControlData","Using control particle set "<<cpii<<" file:"<<cset->mContrPartFile<<" cpmfile:"<<cset->mCpmotionFile<<" mirr:'"<<cset->mCpmotion->getInitMirror()<<"' " );
185                         mCons.push_back( cset );
186                 } else {
187                         delete cset;
188                 }
189         }
190
191         // debug, testing - make sure theres at least an empty set
192         if(mCons.size()<1) {
193                 mCons.push_back( new LbmControlSet() );
194                 mCons[0]->initCparts();
195         }
196
197         // take from first set
198         for(int cpii=1; cpii<(int)mCons.size(); cpii++) {
199                 mCons[cpii]->mCparts->setRadiusMinMaxd( mCons[0]->mCparts->getRadiusMinMaxd() );
200                 mCons[cpii]->mCparts->setRadiusMaxd(    mCons[0]->mCparts->getRadiusMaxd() );
201                 mCons[cpii]->mCparts->setInfluenceAttraction( mCons[0]->mCparts->getInfluenceAttraction() );
202                 mCons[cpii]->mCparts->setInfluenceTangential( mCons[0]->mCparts->getInfluenceTangential() );
203                 mCons[cpii]->mCparts->setInfluenceVelocity( mCons[0]->mCparts->getInfluenceVelocity() , 0.01 ); // dummy dt
204                 mCons[cpii]->mCparts->setInfluenceMaxdist( mCons[0]->mCparts->getInfluenceMaxdist() );
205         }
206         
207         // invert for usage in relax macro
208         mDiffVelCon =  1.-attr->readFloat("cpdiffvelcon",  mDiffVelCon, "LbmControlData","mDiffVelCon", false);
209
210         mDebugCpscale     =  attr->readFloat("cpdebug_cpscale",  mDebugCpscale, "LbmControlData","mDebugCpscale", false);
211         mDebugMaxdScale   =  attr->readFloat("cpdebug_maxdscale",  mDebugMaxdScale, "LbmControlData","mDebugMaxdScale", false);
212         mDebugAttScale    =  attr->readFloat("cpdebug_attscale",  mDebugAttScale, "LbmControlData","mDebugAttScale", false);
213         mDebugVelScale    =  attr->readFloat("cpdebug_velscale",  mDebugVelScale, "LbmControlData","mDebugVelScale", false);
214         mDebugCompavScale =  attr->readFloat("cpdebug_compavscale",  mDebugCompavScale, "LbmControlData","mDebugCompavScale", false);
215         mDebugAvgVelScale =  attr->readFloat("cpdebug_avgvelsc",  mDebugAvgVelScale, "LbmControlData","mDebugAvgVelScale", false);
216 }
217
218
219 void 
220 LbmFsgrSolver::initCpdata()
221 {
222         // enable for cps via env. vars
223         //if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){ mUseTestdata=1; }
224
225         
226         // manually switch on! if this is zero, nothing is done...
227         mpControl->mSetForceStrength = this->mTForceStrength = 1.;
228         mpControl->mCons.clear();
229         
230         // init all control fluid objects
231         int numobjs = (int)(mpGiObjects->size());
232         for(int o=0; o<numobjs; o++) {
233                 ntlGeometryObjModel *obj = (ntlGeometryObjModel *)(*mpGiObjects)[o];
234                 if(obj->getGeoInitType() & FGI_CONTROL) {
235                         // add new control set per object
236                         LbmControlSet *cset;
237
238                         cset = new LbmControlSet();
239                         cset->initCparts();
240         
241                         // dont load any file
242                         cset->mContrPartFile = string("");
243
244                         cset->mcForceAtt = obj->getCpsAttrFStr();
245                         cset->mcRadiusAtt = obj->getCpsAttrFRad();
246                         cset->mcForceVel = obj->getCpsVelFStr();
247                         cset->mcRadiusVel = obj->getCpsVelFRad();
248
249                         cset->mCparts->setCPSTimeStart(obj->getCpsTimeStart());
250                         cset->mCparts->setCPSTimeEnd(obj->getCpsTimeEnd());
251                         
252                         if(obj->getCpsQuality() > LBM_EPSILON)
253                                 cset->mCparts->setCPSWith(1.0 / obj->getCpsQuality());
254                         
255                         // this value can be left at 0.5:
256                         cset->mCparts->setCPSMvmWeightFac(0.5);
257
258                         mpControl->mCons.push_back( cset );
259                         mpControl->mCons[mpControl->mCons.size()-1]->mCparts->initFromObject(obj);
260                 }
261         }
262         
263         // NT blender integration manual test setup
264         if(0) {
265                 // manually switch on! if this is zero, nothing is done...
266                 mpControl->mSetForceStrength = this->mTForceStrength = 1.;
267                 mpControl->mCons.clear();
268
269                 // add new set
270                 LbmControlSet *cset;
271
272                 cset = new LbmControlSet();
273                 cset->initCparts();
274                 // dont load any file
275                 cset->mContrPartFile = string("");
276
277                 // set radii for attraction & velocity forces
278                 // set strength of the forces
279                 // don't set directly! but use channels: 
280                 // mcForceAtt, mcForceVel, mcForceMaxd, mcRadiusAtt, mcRadiusVel, mcRadiusMind, mcRadiusMaxd etc.
281
282                 // wrong: cset->mCparts->setInfluenceAttraction(1.15); cset->mCparts->setRadiusAtt(1.5);
283                 // right, e.g., to init some constant values:
284                 cset->mcForceAtt = AnimChannel<float>(0.2);
285                 cset->mcRadiusAtt = AnimChannel<float>(0.75);
286                 cset->mcForceVel = AnimChannel<float>(0.2);
287                 cset->mcRadiusVel = AnimChannel<float>(0.75);
288
289                 // this value can be left at 0.5:
290                 cset->mCparts->setCPSMvmWeightFac(0.5);
291
292                 mpControl->mCons.push_back( cset );
293
294                 // instead of reading from file (cset->mContrPartFile), manually init some particles
295                 mpControl->mCons[0]->mCparts->initBlenderTest();
296
297                 // other values that might be interesting to change:
298                 //cset->mCparts->setCPSTimestep(0.02);
299                 //cset->mCparts->setCPSTimeStart(0.);
300                 //cset->mCparts->setCPSTimeEnd(1.); 
301
302                 //mpControl->mDiffVelCon = 1.; // more rigid velocity control, 0 (default) allows more turbulence
303         }
304
305         // control particle -------------------------------------------------------------------------------------
306
307         // init cppf stage, use set 0!
308         if(mCppfStage>0) {
309                 if(mpControl->mCpOutfile.length()<1) mpControl->mCpOutfile = string("cpout"); // use getOutFilename !?
310                 char strbuf[100];
311                 const char *cpFormat = "_d%dcppf%d";
312
313                 // initial coarse stage, no input
314                 if(mCppfStage==1) {
315                         mpControl->mCons[0]->mContrPartFile = "";
316                 } else {
317                         // read from prev stage
318                         snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage-1);
319                         mpControl->mCons[0]->mContrPartFile = mpControl->mCpOutfile;
320                         mpControl->mCons[0]->mContrPartFile += strbuf;
321                         mpControl->mCons[0]->mContrPartFile += ".cpart2";
322                 }
323
324                 snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage);
325                 mpControl->mCpOutfile += strbuf;
326         } // */
327         
328         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
329                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
330                 ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
331
332                 // now set with real dt
333                 cparts->setInfluenceVelocity( mpControl->mCons[cpssi]->mcForceVel.get(0.), mLevel[mMaxRefine].timestep);
334                 cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
335                 cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
336                 errMsg("LbmControlData","CppfStage "<<mCppfStage<<" in:"<<mpControl->mCons[cpssi]->mContrPartFile<<
337                                 " out:"<<mpControl->mCpOutfile<<" cl:"<< cparts->getCharLength() );
338
339                 // control particle test init
340                 if(mpControl->mCons[cpssi]->mCpmotionFile.length()>=1) cpmotion->initFromTextFile(mpControl->mCons[cpssi]->mCpmotionFile);
341                 // not really necessary...
342                 //? cparts->setFluidSpacing( mLevel[mMaxRefine].nodeSize        ); // use grid coords!?
343                 //? cparts->calculateKernelWeight();
344                 //? debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - motion inited: "<<cparts->getSize() ,10);
345
346                 // ensure both are on for env. var settings
347                 // when no particles, but outfile enabled, initialize
348                 const int lev = mMaxRefine;
349                 if((mpParticles) && (mpControl->mCpOutfile.length()>=1) && (cpssi==0)) {
350                         // check if auto num
351                         if( (mpParticles->getNumInitialParticles()<=1) && 
352                                         (mpParticles->getNumParticles()<=1) ) { // initParticles done afterwards anyway
353                                 int tracers = 0;
354                                 const int workSet = mLevel[lev].setCurr;
355                                 FSGR_FORIJK_BOUNDS(lev) { 
356                                         if(RFLAG(lev,i,j,k, workSet)&(CFFluid)) tracers++;
357                                 }
358                                 if(LBMDIM==3) tracers /= 8;
359                                 else          tracers /= 4;
360                                 mpParticles->setNumInitialParticles(tracers);
361                                 mpParticles->setDumpTextFile(mpControl->mCpOutfile);
362                                 debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - set tracers #"<<tracers<<", actual #"<<mpParticles->getNumParticles() ,10);
363                         }
364                         if(mpParticles->getDumpTextInterval()<=0.) {
365                                 mpParticles->setDumpTextInterval(mLevel[lev].timestep * mLevel[lev].lSizex);
366                                 debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - dump delta t not set, using dti="<< mpParticles->getDumpTextInterval()<<", sim dt="<<mLevel[lev].timestep, 5 );
367                         }
368                         mpParticles->setDumpParts(true); // DEBUG? also dump as particle system
369                 }
370
371                 if(mpControl->mCons[cpssi]->mContrPartFile.length()>=1) cparts->initFromTextFile(mpControl->mCons[cpssi]->mContrPartFile);
372                 cparts->setFluidSpacing( mLevel[lev].nodeSize   ); // use grid coords!?
373                 cparts->calculateKernelWeight();
374                 debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles mCons"<<cpssi<<" - inited, parts:"<<cparts->getTotalSize()<<","<<cparts->getSize()<<" dt:"<<mpParam->getTimestep()<<" control time:"<<cparts->getControlTimStart()<<" to "<<cparts->getControlTimEnd() ,10);
375         } // cpssi
376
377         if(getenv("ELBEEM_CPINFILE")) {
378                 this->mTForceStrength = 1.0;
379         }
380         this->mTForceStrength = mpControl->mSetForceStrength;
381         if(mpControl->mCpOutfile.length()>=1) mpParticles->setDumpTextFile(mpControl->mCpOutfile);
382
383         // control particle  init end -------------------------------------------------------------------------------------
384
385         // make sure equiv to solver init
386         if(this->mTForceStrength>0.) { \
387                 mpControl->mCpForces.resize( mMaxRefine+1 );
388                 for(int lev = 0; lev<=mMaxRefine; lev++) {
389                         LONGINT rcellSize = (mLevel[lev].lSizex*mLevel[lev].lSizey*mLevel[lev].lSizez);
390                         debMsgStd("LbmFsgrSolver::initControl",DM_MSG,"mCpForces init, lev="<<lev<<" rcs:"<<(int)(rcellSize+4)<<","<<(rcellSize*sizeof(ControlForces)/(1024*1024)), 9 );
391                         mpControl->mCpForces[lev].resize( (int)(rcellSize+4) );
392                         //for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces.push_back( ControlForces() );
393                         for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces[lev][i].resetForces();
394                 }
395         } // on?
396
397         debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles #mCons "<<mpControl->mCons.size()<<" done", 6);
398 }
399
400
401 #define CPODEBUG 0
402 //define CPINTER ((int)(mpControl->mCpUpdateInterval))
403
404 #define KERN(x,y,z)    mpControl->mCpKernel[ (((z)*cpkarWidth + (y))*cpkarWidth + (x)) ]
405 #define MDKERN(x,y,z)  mpControl->mMdKernel[ (((z)*mdkarWidth + (y))*mdkarWidth + (x)) ]
406
407 #define BOUNDCHECK(x,low,high)  ( ((x)<low) ? low : (((x)>high) ? high : (x) )  )
408 #define BOUNDSKIP(x,low,high)  ( ((x)<low) || ((x)>high) )
409
410 void 
411 LbmFsgrSolver::handleCpdata()
412 {
413         myTime_t cpstart = getTime();
414         int cpChecks=0;
415         int cpInfs=0;
416         //debMsgStd("ControlData::handleCpdata",DM_MSG,"called... "<<this->mTForceStrength,1);
417
418         // add cp influence
419   if((true) && (this->mTForceStrength>0.)) {
420                 // ok continue...
421         } // on off
422         else {
423                 return;
424         }
425         
426         // check if we have control objects
427         if(mpControl->mCons.size()==0)
428                 return;
429         
430         if((mpControl->mCpUpdateInterval<1) || (this->mStepCnt%mpControl->mCpUpdateInterval==0)) {
431                 // do full reinit later on... 
432         }
433         else if(this->mStepCnt>mpControl->mCpUpdateInterval) {
434                 // only reinit new cells
435                 // TODO !? remove loop dependance!?
436 #define NOFORCEENTRY(lev, i,j,k) (LBMGET_FORCE(lev, i,j,k).maxDistance==CPF_MAXDINIT) 
437                 // interpolate missing
438                 for(int lev=0; lev<=mMaxRefine; lev++) {
439                 FSGR_FORIJK_BOUNDS(lev) { 
440                         if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFFluid|CFInter) )  
441                         //if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFInter) )  
442                         //if(0)
443                         { // only check new inter? RFLAG?check
444                                 if(NOFORCEENTRY(lev, i,j,k)) {
445                                         //errMsg("CP","FE_MISSING at "<<PRINT_IJK<<" f"<<LBMGET_FORCE(lev, i,j,k).weightAtt<<" md"<<LBMGET_FORCE(lev, i,j,k).maxDistance );
446
447                                         LbmFloat nbs=0.;
448                                         ControlForces vals;
449                                         vals.resetForces(); vals.maxDistance = 0.;
450                                         for(int l=1; l<this->cDirNum; l++) { 
451                                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
452                                                 //errMsg("CP","FE_MISSING check "<<PRINT_VEC(ni,nj,nk)<<" f"<<LBMGET_FORCE(lev, ni,nj,nk).weightAtt<<" md"<<LBMGET_FORCE(lev, ni,nj,nk).maxDistance );
453                                                 if(!NOFORCEENTRY(lev, ni,nj,nk)) {
454                                                         //? vals.weightAtt   += LBMGET_FORCE(lev, ni,nj,nk).weightAtt;
455                                                         //? vals.forceAtt    += LBMGET_FORCE(lev, ni,nj,nk).forceAtt;
456                                                         vals.maxDistance += LBMGET_FORCE(lev, ni,nj,nk).maxDistance;
457                                                         vals.forceMaxd   += LBMGET_FORCE(lev, ni,nj,nk).forceMaxd;
458                                                         vals.weightVel   += LBMGET_FORCE(lev, ni,nj,nk).weightVel;
459                                                         vals.forceVel    += LBMGET_FORCE(lev, ni,nj,nk).forceVel; 
460                                                         // ignore att/compAv/avgVel here for now
461                                                         nbs += 1.;
462                                                 }
463                                         }
464                                         if(nbs>0.) {
465                                                 nbs = 1./nbs;
466                                                 //? LBMGET_FORCE(lev, i,j,k).weightAtt   =      vals.weightAtt*nbs;
467                                                 //? LBMGET_FORCE(lev, i,j,k).forceAtt    = vals.forceAtt*nbs;
468                                                 LBMGET_FORCE(lev, i,j,k).maxDistance = vals.maxDistance*nbs;
469                                                 LBMGET_FORCE(lev, i,j,k).forceMaxd   =  vals.forceMaxd*nbs;
470                                                 LBMGET_FORCE(lev, i,j,k).weightVel   =  vals.weightVel*nbs;
471                                                 LBMGET_FORCE(lev, i,j,k).forceVel    = vals.forceVel*nbs;
472                                         }
473                                                 /*ControlForces *ff = &LBMGET_FORCE(lev, i,j,k);  // DEBUG
474                                                 errMsg("CP","FE_MISSING rec at "<<PRINT_IJK // DEBUG
475                                                                 <<" w:"<<ff->weightAtt<<" wa:" <<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
476                                                                 <<" v:"<<ff->weightVel<<" wv:" <<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
477                                                                 <<" v:"<<ff->maxDistance<<" wv:" <<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  ); // DEBUG */
478                                         // else errMsg("CP","FE_MISSING rec at "<<PRINT_IJK<<" failed!"); // DEBUG
479                                         
480                                 }
481                         }
482                 }} // ijk, lev
483
484                 // mStepCnt > mpControl->mCpUpdateInterval
485                 return;
486         } else {
487                 // nothing to do ...
488                 return;
489         }
490
491         // reset
492         for(int lev=0; lev<=mMaxRefine; lev++) {
493                 FSGR_FORIJK_BOUNDS(lev) { LBMGET_FORCE(lev,i,j,k).resetForces(); }
494         }
495         // do setup for coarsest level
496         const int coarseLev = 0;
497         const int fineLev = mMaxRefine;
498
499         // init for current time
500         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
501                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
502                 LbmControlSet *cset = mpControl->mCons[cpssi];
503                 
504                 cparts->setRadiusAtt(cset->mcRadiusAtt.get(mSimulationTime));
505                 cparts->setRadiusVel(cset->mcRadiusVel.get(mSimulationTime));
506                 cparts->setInfluenceAttraction(cset->mcForceAtt.get(mSimulationTime) );
507                 cparts->setInfluenceMaxdist(cset->mcForceMaxd.get(mSimulationTime) );
508                 cparts->setRadiusMinMaxd(cset->mcRadiusMind.get(mSimulationTime));
509                 cparts->setRadiusMaxd(cset->mcRadiusMaxd.get(mSimulationTime));
510                 cparts->calculateKernelWeight(); // always necessary!?
511                 cparts->setOffset( vec2L(cset->mcCpOffset.get(mSimulationTime)) );
512                 cparts->setScale(  vec2L(cset->mcCpScale.get(mSimulationTime)) );
513
514                 cparts->setInfluenceVelocity( cset->mcForceVel.get(mSimulationTime), mLevel[fineLev].timestep );
515                 cparts->setLastOffset( vec2L(cset->mcCpOffset.get(mSimulationTime-mLevel[fineLev].timestep)) );
516                 cparts->setLastScale(  vec2L(cset->mcCpScale.get(mSimulationTime-mLevel[fineLev].timestep)) );
517                 
518         }
519
520         // check actual values
521         LbmFloat iatt  = ABS(mpControl->mCons[0]->mCparts->getInfluenceAttraction());
522         LbmFloat ivel  = mpControl->mCons[0]->mCparts->getInfluenceVelocity();
523         LbmFloat imaxd = mpControl->mCons[0]->mCparts->getInfluenceMaxdist();
524         //errMsg("FINCIT","iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
525         for(int cpssi=1; cpssi<(int)mpControl->mCons.size(); cpssi++) {
526                 LbmFloat iatt2  = ABS(mpControl->mCons[cpssi]->mCparts->getInfluenceAttraction());
527                 LbmFloat ivel2  = mpControl->mCons[cpssi]->mCparts->getInfluenceVelocity();
528                 LbmFloat imaxd2 = mpControl->mCons[cpssi]->mCparts->getInfluenceMaxdist();
529                 
530                 // we allow negative attraction force here!
531                 if(iatt2 > iatt)   iatt = iatt2;
532                 
533                 if(ivel2 >ivel)   ivel = ivel2;
534                 if(imaxd2>imaxd)  imaxd= imaxd2;
535                 //errMsg("FINCIT"," "<<cpssi<<" iatt2="<<iatt2<<" ivel2="<<ivel2<<" imaxd2="<<imaxd<<" NEW "<<" iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
536         }
537
538         if(iatt==0. && ivel==0. && imaxd==0.) {
539                 debMsgStd("ControlData::initControl",DM_MSG,"Skipped, all zero...",4);
540                 return;
541         }
542         //iatt  = mpControl->mCons[1]->mCparts->getInfluenceAttraction(); //ivel  = mpControl->mCons[1]->mCparts->getInfluenceVelocity(); //imaxd = mpControl->mCons[1]->mCparts->getInfluenceMaxdist(); // TTTTTT
543
544         // do control setup
545         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
546                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
547                 ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
548
549                 // TEST!?
550                 bool radmod = false;
551                 const LbmFloat minRadSize = mLevel[coarseLev].nodeSize * 1.5;
552                 if((cparts->getRadiusAtt()>0.) && (cparts->getRadiusAtt()<minRadSize) && (!radmod) ) {
553                         LbmFloat radfac = minRadSize / cparts->getRadiusAtt(); radmod=true;
554                         debMsgStd("ControlData::initControl",DM_MSG,"Modified radii att, fac="<<radfac, 7);
555                         cparts->setRadiusAtt(cparts->getRadiusAtt()*radfac);
556                         cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
557                         cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
558                         cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
559                 } else if((cparts->getRadiusVel()>0.) && (cparts->getRadiusVel()<minRadSize) && (!radmod) ) {
560                         LbmFloat radfac = minRadSize / cparts->getRadiusVel();
561                         debMsgStd("ControlData::initControl",DM_MSG,"Modified radii vel, fac="<<radfac, 7);
562                         cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
563                         cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
564                         cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
565                 } else if((cparts->getRadiusMaxd()>0.) && (cparts->getRadiusMaxd()<minRadSize) && (!radmod) ) {
566                         LbmFloat radfac = minRadSize / cparts->getRadiusMaxd();
567                         debMsgStd("ControlData::initControl",DM_MSG,"Modified radii maxd, fac="<<radfac, 7);
568                         cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
569                         cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
570                 }
571                 if(radmod) {
572                         debMsgStd("ControlData::initControl",DM_MSG,"Modified radii: att="<< 
573                                 cparts->getRadiusAtt()<<", vel=" << cparts->getRadiusVel()<<", maxd=" <<
574                                 cparts->getRadiusMaxd()<<", mind=" << cparts->getRadiusMinMaxd() ,5);
575                 }
576
577                 cpmotion->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), NULL );
578                 cparts->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), cpmotion );
579         }
580
581         // do control...
582         for(int lev=0; lev<=mMaxRefine; lev++) {
583                 LbmFloat levVolume = 1.;
584                 LbmFloat levForceScale = 1.;
585                 for(int ll=lev; ll<mMaxRefine; ll++) {
586                         if(LBMDIM==3) levVolume *= 8.;
587                         else levVolume *= 4.;
588                         levForceScale *= 2.;
589                 }
590                 errMsg("LbmFsgrSolver::handleCpdata","levVolume="<<levVolume<<" levForceScale="<<levForceScale );
591         //todo: scale velocity, att by level timestep!?
592                 
593         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
594                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
595                 // ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
596                 
597                 // if control set is not active skip it
598                 if((cparts->getControlTimStart() > mSimulationTime) || (cparts->getControlTimEnd() < mLastSimTime))
599                 {
600                         continue;
601                 }
602
603                 const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
604                 LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
605                 LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
606                 LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
607 #if LBMDIM==2
608                 gsz = gsx;
609 #endif
610                 LbmFloat goffx = mvGeoStart[0];
611                 LbmFloat goffy = mvGeoStart[1];
612                 LbmFloat goffz = mvGeoStart[2];
613
614                 //const LbmFloat cpwIncFac = 2.0;
615                 // max to two thirds of domain size
616                 const int cpw = MIN( mLevel[lev].lSizex/3, MAX( (int)( cparts->getRadiusAtt()  /gsx) +1  , 2) ); // normal kernel, att,vel
617                 const int cpkarWidth = 2*cpw+1;
618                 mpControl->mCpKernel.resize(cpkarWidth* cpkarWidth* cpkarWidth);
619                 ControlParticle cpt; cpt.reset();
620                 cpt.pos = LbmVec( (gsx*(LbmFloat)cpw)+goffx, (gsy*(LbmFloat)cpw)+goffy, (gsz*(LbmFloat)cpw)+goffz );  // optimize?
621                 cpt.density = 0.5; cpt.densityWeight = 0.5;
622 #if LBMDIM==3
623                 for(int k= 0; k<cpkarWidth; ++k) {
624 #else // LBMDIM==3
625                 { int k = cpw;
626 #endif 
627                         for(int j= 0; j<cpkarWidth; ++j) 
628                                 for(int i= 0; i<cpkarWidth; ++i) {
629                                         KERN(i,j,k).resetForces();
630                                         //LbmFloat dx = i-cpw; LbmFloat dy = j-cpw; LbmFloat dz = k-cpw;
631                                         //LbmVec dv = ( LbmVec(dx,dy,dz) );
632                                         //LbmFloat dl = norm( dv ); //LbmVec dir = dv / dl;
633                                         LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
634                                         cparts->calculateCpInfluenceOpt( &cpt, pos, LbmVec(0,0,0), &KERN(i,j,k)  ,1. );
635                                         /*if((CPODEBUG)&&(k==cpw)) errMsg("kern"," at "<<PRINT_IJK<<" pos"<<pos<<" cpp"<<cpt.pos 
636                                                 <<" wf:"<<KERN(i,j,k).weightAtt<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceAtt[0],KERN(i,j,k).forceAtt[1],KERN(i,j,k).forceAtt[2] )  
637                                                 <<" wf:"<<KERN(i,j,k).weightVel<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceVel[0],KERN(i,j,k).forceVel[1],KERN(i,j,k).forceVel[2] )  
638                                                 <<" wf:"<<KERN(i,j,k).maxDistance<<" wa:"<< PRINT_VEC( KERN(i,j,k).forceMaxd[0],KERN(i,j,k).forceMaxd[1],KERN(i,j,k).forceMaxd[2] )  ); // */
639                                         KERN(i,j,k).weightAtt *= 2.;
640                                         KERN(i,j,k).forceAtt *= 2.;
641                                         //KERN(i,j,k).forceAtt[1] *= 2.; KERN(i,j,k).forceAtt[2] *= 2.;
642                                         KERN(i,j,k).weightVel *= 2.;
643                                         KERN(i,j,k).forceVel *= 2.;
644                                         //KERN(i,j,k).forceVel[1] *= 2.; KERN(i,j,k).forceVel[2] *= 2.;
645                                 }
646                 }
647
648                 if(CPODEBUG) errMsg("cpw"," = "<<cpw<<" f"<< cparts->getRadiusAtt()<<" gsx"<<gsx<<" kpw"<<cpkarWidth); // DEBUG
649                 // first cp loop - add att and vel forces
650                 for(int cppi=0; cppi<cparts->getSize(); cppi++) {
651                         ControlParticle *cp = cparts->getParticle(cppi);
652                         if(cp->influence<=0.) continue;
653                         const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
654                         const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
655                         int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
656                         /*if( ((LBMDIM==3)&&(BOUNDSKIP(cpk - cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
657                                 ((LBMDIM==3)&&(BOUNDSKIP(cpk + cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
658                                 BOUNDSKIP(cpj - cpwsm, 0, mLevel[lev].lSizey ) ||
659                                 BOUNDSKIP(cpj + cpwsm, 0, mLevel[lev].lSizey ) ||
660                                 BOUNDSKIP(cpi - cpwsm, 0, mLevel[lev].lSizex ) ||
661                                 BOUNDSKIP(cpi + cpwsm, 0, mLevel[lev].lSizex ) ) {
662                                 continue;
663                                 } // */
664                         int is,ie,js,je,ks,ke;
665                         ks = BOUNDCHECK(cpk - cpw, getForZMinBnd(), getForZMaxBnd(lev) );
666                         ke = BOUNDCHECK(cpk + cpw, getForZMinBnd(), getForZMaxBnd(lev) );
667                         js = BOUNDCHECK(cpj - cpw, 0, mLevel[lev].lSizey );
668                         je = BOUNDCHECK(cpj + cpw, 0, mLevel[lev].lSizey );
669                         is = BOUNDCHECK(cpi - cpw, 0, mLevel[lev].lSizex );
670                         ie = BOUNDCHECK(cpi + cpw, 0, mLevel[lev].lSizex );
671                         if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
672                         if(CPODEBUG) errMsg("cppft","i"<<cppi<<" cpw"<<cpw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
673                         cpInfs++;
674
675                         for(int k= ks; k<ke; ++k) {
676                                 for(int j= js; j<je; ++j) {
677
678                                         CellFlagType *pflag = &RFLAG(lev,is,j,k, mLevel[lev].setCurr);
679                                         ControlForces *kk = &KERN( is-cpi+cpw, j-cpj+cpw, k-cpk+cpw);
680                                         ControlForces *ff = &LBMGET_FORCE(lev,is,j,k); 
681                                         pflag--; kk--; ff--;
682
683                                         for(int i= is; i<ie; ++i) {
684                                                 // first cp loop (att,vel)
685                                                 pflag++; kk++; ff++;
686
687                                                 //add weight for bnd cells
688                                                 const LbmFloat pwforce = kk->weightAtt;
689                                                 // control particle mod,
690                                                 // dont add multiple CFFluid fsgr boundaries
691                                                 if(lev==mMaxRefine) {
692                                                         //if( ( ((*pflag)&(CFFluid )) && (lev==mMaxRefine) ) ||
693                                                                         //( ((*pflag)&(CFGrNorm)) && (lev <mMaxRefine) ) ) {
694                                                         if((*pflag)&(CFFluid|CFUnused)) {
695                                                                 // check not fromcoarse?
696                                                                 cp->density += levVolume* kk->weightAtt; // old CFFluid
697                                                         } else if( (*pflag) & (CFEmpty) ) {  
698                                                                 cp->density -= levVolume* 0.5; 
699                                                         } else { //if( ((*pflag) & (CFBnd)) ) {  
700                                                                 cp->density -= levVolume* 0.2;  // penalty
701                                                         }
702                                                 } else {
703                                                         //if((*pflag)&(CFGrNorm)) {
704                                                                 //cp->density += levVolume* kk->weightAtt; // old CFFluid
705                                                         //} 
706                                                 }
707                                                 //else if(!((*pflag) & (CFUnused)) ) {  cp->density -= levVolume* 0.2; } // penalty
708
709                                                 if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
710                                                 {
711
712                                                         cpChecks++;
713                                                         //const LbmFloat pwforce = kk->weightAtt;
714                                                         LbmFloat pwvel = kk->weightVel;
715                                                         if((pwforce==0.)&&(pwvel==0.)) { continue; }
716                                                         ff->weightAtt += 1e-6; // for distance
717
718                                                         if(pwforce>0.) {
719                                                                 ff->weightAtt += pwforce *cp->densityWeight *cp->influence;
720                                                                 ff->forceAtt += kk->forceAtt *levForceScale *cp->densityWeight *cp->influence;
721
722                                                                 // old fill handling here
723                                                                 const int workSet =mLevel[lev].setCurr;
724                                                                 LbmFloat ux=0., uy=0., uz=0.;
725                                                                 FORDF1{  
726                                                                         const LbmFloat dfn = QCELL(lev, i,j,k, workSet, l);
727                                                                         ux  += (this->dfDvecX[l]*dfn);
728                                                                         uy  += (this->dfDvecY[l]*dfn); 
729                                                                         uz  += (this->dfDvecZ[l]*dfn); 
730                                                                 }
731                                                                 // control particle mod
732                                                                 cp->avgVelWeight += levVolume*pwforce;
733                                                                 cp->avgVelAcc += LbmVec(ux,uy,uz) * levVolume*pwforce;
734                                                         }
735
736                                                         if(pwvel>0.) {
737                                                                 // TODO make switch? vel.influence depends on density weight... 
738                                                                 // (reduced lowering with 0.75 factor)
739                                                                 pwvel *=  cp->influence *(1.-0.75*cp->densityWeight);
740                                                                 // control particle mod
741                                                                 // todo use Omega instead!?
742                                                                 ff->forceVel += cp->vel*levVolume*pwvel * velLatticeScale; // levVolume?
743                                                                 ff->weightVel += levVolume*pwvel; // levVolume?
744                                                                 ff->compAv += cp->avgVel*levVolume*pwvel; // levVolume?
745                                                                 ff->compAvWeight += levVolume*pwvel; // levVolume?
746                                                         }
747
748                                                         if(CPODEBUG) errMsg("cppft","i"<<cppi<<" at "<<PRINT_IJK<<" kern:"<<
749                                                                         PRINT_VEC(i-cpi+cpw, j-cpj+cpw, k-cpk+cpw )
750                                                                         //<<" w:"<<ff->weightAtt<<" wa:"
751                                                                         //<<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
752                                                                         //<<" v:"<<ff->weightVel<<" wv:"
753                                                                         //<<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
754                                                                         //<<" v:"<<ff->maxDistance<<" wv:"
755                                                                         //<<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  
756                                                                         );
757                                                 } // celltype
758
759                                         } // ijk
760                                 } // ijk
761                         } // ijk
762                 } // cpi, end first cp loop (att,vel)
763                 debMsgStd("LbmFsgrSolver::handleCpdata",DM_MSG,"Force cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
764         } //cpssi
765         } // lev
766
767         // second loop
768         for(int lev=0; lev<=mMaxRefine; lev++) {
769                 LbmFloat levVolume = 1.;
770                 LbmFloat levForceScale = 1.;
771                 for(int ll=lev; ll<mMaxRefine; ll++) {
772                         if(LBMDIM==3) levVolume *= 8.;
773                         else levVolume *= 4.;
774                         levForceScale *= 2.;
775                 }
776         // prepare maxd forces
777         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
778                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
779
780                 // WARNING copied from above!
781                 const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
782                 LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
783                 LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
784                 LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
785 #if LBMDIM==2
786                 gsz = gsx;
787 #endif
788                 LbmFloat goffx = mvGeoStart[0];
789                 LbmFloat goffy = mvGeoStart[1];
790                 LbmFloat goffz = mvGeoStart[2];
791
792                 //const LbmFloat cpwIncFac = 2.0;
793                 const int mdw = MIN( mLevel[lev].lSizex/2, MAX( (int)( cparts->getRadiusMaxd() /gsx) +1  , 2) ); // wide kernel, md
794                 const int mdkarWidth = 2*mdw+1;
795                 mpControl->mMdKernel.resize(mdkarWidth* mdkarWidth* mdkarWidth);
796                 ControlParticle cpt; cpt.reset();
797                 cpt.density = 0.5; cpt.densityWeight = 0.5;
798                 cpt.pos = LbmVec( (gsx*(LbmFloat)mdw)+goffx, (gsy*(LbmFloat)mdw)+goffy, (gsz*(LbmFloat)mdw)+goffz );  // optimize?
799 #if LBMDIM==3
800                 for(int k= 0; k<mdkarWidth; ++k) {
801 #else // LBMDIM==3
802                 { int k = mdw;
803 #endif 
804                         for(int j= 0; j<mdkarWidth; ++j) 
805                                 for(int i= 0; i<mdkarWidth; ++i) {
806                                         MDKERN(i,j,k).resetForces();
807                                         LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
808                                         cparts->calculateMaxdForce( &cpt, pos, &MDKERN(i,j,k)  );
809                                 }
810                 }
811
812                 // second cpi loop, maxd forces
813                 if(cparts->getInfluenceMaxdist()>0.) {
814                         for(int cppi=0; cppi<cparts->getSize(); cppi++) {
815                                 ControlParticle *cp = cparts->getParticle(cppi);
816                                 if(cp->influence<=0.) continue;
817                                 const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
818                                 const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
819                                 int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
820
821                                 int is,ie,js,je,ks,ke;
822                                 ks = BOUNDCHECK(cpk - mdw, getForZMinBnd(), getForZMaxBnd(lev) );
823                                 ke = BOUNDCHECK(cpk + mdw, getForZMinBnd(), getForZMaxBnd(lev) );
824                                 js = BOUNDCHECK(cpj - mdw, 0, mLevel[lev].lSizey );
825                                 je = BOUNDCHECK(cpj + mdw, 0, mLevel[lev].lSizey );
826                                 is = BOUNDCHECK(cpi - mdw, 0, mLevel[lev].lSizex );
827                                 ie = BOUNDCHECK(cpi + mdw, 0, mLevel[lev].lSizex );
828                                 if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
829                                 if(CPODEBUG) errMsg("cppft","i"<<cppi<<" mdw"<<mdw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
830                                 cpInfs++;
831
832                                 for(int k= ks; k<ke; ++k)
833                                  for(int j= js; j<je; ++j) {
834                                         CellFlagType *pflag = &RFLAG(lev,is-1,j,k, mLevel[lev].setCurr);
835                                         for(int i= is; i<ie; ++i) {
836                                                 // second cpi loop, maxd forces
837                                                 pflag++;
838                                                 if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
839                                                 {
840                                                         cpChecks++;
841                                                         ControlForces *ff = &LBMGET_FORCE(lev,i,j,k); 
842                                                         if(ff->weightAtt == 0.) {
843                                                                 ControlForces *kk = &MDKERN( i-cpi+mdw, j-cpj+mdw, k-cpk+mdw);
844                                                                 const LbmFloat pmdf = kk->maxDistance;
845                                                                 if((ff->maxDistance > pmdf) || (ff->maxDistance<0.))
846                                                                         ff->maxDistance = pmdf;
847                                                                 ff->forceMaxd = kk->forceMaxd;
848                                                                 // todo use Omega instead!?
849                                                                 ff->forceVel = cp->vel* velLatticeScale;
850                                                         }
851                                                 } // celltype
852                                 } } // ijk
853                         } // cpi, md loop 
854                 } // maxd inf>0 */
855
856
857                 debMsgStd("ControlData::initControl",DM_MSG,"Maxd cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
858         } //cpssi
859
860         // normalize, only done once for the whole array
861         mpControl->mCons[0]->mCparts->finishControl( mpControl->mCpForces[lev], iatt,ivel,imaxd );
862
863         } // lev loop
864
865         myTime_t cpend = getTime(); 
866         debMsgStd("ControlData::handleCpdata",DM_MSG,"Time for cpgrid generation:"<< getTimeString(cpend-cpstart)<<", checks:"<<cpChecks<<" infs:"<<cpInfs<<" " ,8);
867
868         // warning, may return  before
869 }
870
871 #if LBM_USE_GUI==1
872
873 #define USE_GLUTILITIES
874 #include "../gui/gui_utilities.h"
875
876 void LbmFsgrSolver::cpDebugDisplay(int dispset)
877 {
878         for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
879                 ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
880                 //ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
881                 // display cp parts
882                 const bool cpCubes = false;
883                 const bool cpDots = true;
884                 const bool cpCpdist = true;
885                 const bool cpHideIna = true;
886                 glShadeModel(GL_FLAT);
887                 glDisable( GL_LIGHTING ); // dont light lines
888
889                 // dot influence
890                 if((mpControl->mDebugCpscale>0.) && cpDots) {
891                         glPointSize(mpControl->mDebugCpscale * 8.);
892                         glBegin(GL_POINTS);
893                         for(int i=0; i<cparts->getSize(); i++) {
894                                 if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
895                                 ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
896                                 //LbmFloat halfsize = 0.5; 
897                                 LbmFloat scale = cparts->getParticle(i)->densityWeight;
898                                 //glColor4f( scale,scale,scale,scale );
899                                 glColor4f( 0.,scale,0.,scale );
900                                 glVertex3f( org[0],org[1],org[2] );
901                                 //errMsg("lbmDebugDisplay","CP "<<i<<" at "<<org); // DEBUG
902                         }
903                         glEnd();
904                 }
905
906                 // cp positions
907                 if((mpControl->mDebugCpscale>0.) && cpDots) {
908                         glPointSize(mpControl->mDebugCpscale * 3.);
909                         glBegin(GL_POINTS); 
910                         glColor3f( 0,1,0 );
911                 }
912                 for(int i=0; i<cparts->getSize(); i++) {
913                         if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
914                         ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
915                         LbmFloat halfsize = 0.5; 
916                         LbmFloat scale = cparts->getRadiusAtt() * cparts->getParticle(i)->densityWeight;
917                         if(cpCubes){    glLineWidth( 1 ); 
918                                 glColor3f( 1,1,1 );
919                                 ntlVec3Gfx s = org-(halfsize * (scale)); 
920                                 ntlVec3Gfx e = org+(halfsize * (scale)); 
921                                 drawCubeWire( s,e ); }
922                         if((mpControl->mDebugCpscale>0.) && cpDots) {
923                                 glVertex3f( org[0],org[1],org[2] );
924                         }
925                 }
926                 if(cpDots) glEnd();
927
928                 if(mpControl->mDebugAvgVelScale>0.) {
929                         const float scale = mpControl->mDebugAvgVelScale;
930
931                         glColor3f( 1.0,1.0,1 );
932                         glBegin(GL_LINES); 
933                         for(int i=0; i<cparts->getSize(); i++) {
934                                 if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
935                                 ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
936
937                                 //errMsg("CPAVGVEL","i"<<i<<" pos"<<org<<" av"<<cparts->getParticle(i)->avgVel);// DEBUG
938                                 float dx = cparts->getParticle(i)->avgVel[0];
939                                 float dy = cparts->getParticle(i)->avgVel[1];
940                                 float dz = cparts->getParticle(i)->avgVel[2];
941                                 dx *= scale; dy *= scale; dz *= scale;
942                                 glVertex3f( org[0],org[1],org[2] );
943                                 glVertex3f( org[0]+dx,org[1]+dy,org[2]+dz );
944                         }
945                         glEnd();
946                 } // */
947
948                 if( (LBMDIM==2) && (cpCpdist) ) {
949                         
950                         // debug, for use of e.g. LBMGET_FORCE LbmControlData *mpControl = this;
951 #                       define TESTGET_FORCE(lev,i,j,k)   mpControl->mCpForces[lev][ ((k*mLevel[lev].lSizey)+j)*mLevel[lev].lSizex+i ]
952         
953                         glBegin(GL_LINES);
954                         //const int lev=0; 
955                         for(int lev=0; lev<=mMaxRefine; lev++) {
956                         FSGR_FORIJK_BOUNDS(lev) {
957                                         LbmVec pos = LbmVec( 
958                                                 ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex) * ((LbmFloat)i+0.5) + mvGeoStart[0], 
959                                                 ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey) * ((LbmFloat)j+0.5) + mvGeoStart[1], 
960                                                 ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez) * ((LbmFloat)k+0.5) + mvGeoStart[2]  );
961                                         if(LBMDIM==2) pos[2] = ((mvGeoEnd[2]-mvGeoStart[2])*0.5 + mvGeoStart[2]);
962
963                                         if((mpControl->mDebugMaxdScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt<=0.) )
964                                         if(TESTGET_FORCE(lev,i,j,k).maxDistance>=0.) 
965                                         if(TESTGET_FORCE(lev,i,j,k).maxDistance<CPF_MAXDINIT ) {
966                                                 const float scale = mpControl->mDebugMaxdScale*10001.;
967                                                 float dx = TESTGET_FORCE(lev,i,j,k).forceMaxd[0];
968                                                 float dy = TESTGET_FORCE(lev,i,j,k).forceMaxd[1];
969                                                 float dz = TESTGET_FORCE(lev,i,j,k).forceMaxd[2];
970                                                 dx *= scale; dy *= scale; dz *= scale;
971                                                 glColor3f( 0,1,0 );
972                                                 glVertex3f( pos[0],pos[1],pos[2] );
973                                                 glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
974                                         } // */
975                                         if((mpControl->mDebugAttScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt>0.)) {
976                                                 const float scale = mpControl->mDebugAttScale*100011.;
977                                                 float dx = TESTGET_FORCE(lev,i,j,k).forceAtt[0];
978                                                 float dy = TESTGET_FORCE(lev,i,j,k).forceAtt[1];
979                                                 float dz = TESTGET_FORCE(lev,i,j,k).forceAtt[2];
980                                                 dx *= scale; dy *= scale; dz *= scale;
981                                                 glColor3f( 1,0,0 );
982                                                 glVertex3f( pos[0],pos[1],pos[2] );
983                                                 glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
984                                         } // */
985                                         // why check maxDistance?
986                                         if((mpControl->mDebugVelScale>0.) && (TESTGET_FORCE(lev,i,j,k).maxDistance+TESTGET_FORCE(lev,i,j,k).weightVel>0.)) {
987                                                 float scale = mpControl->mDebugVelScale*1.;
988                                                 float wvscale = TESTGET_FORCE(lev,i,j,k).weightVel;
989                                                 float dx = TESTGET_FORCE(lev,i,j,k).forceVel[0];
990                                                 float dy = TESTGET_FORCE(lev,i,j,k).forceVel[1];
991                                                 float dz = TESTGET_FORCE(lev,i,j,k).forceVel[2];
992                                                 scale *= wvscale;
993                                                 dx *= scale; dy *= scale; dz *= scale;
994                                                 glColor3f( 0.2,0.2,1 );
995                                                 glVertex3f( pos[0],pos[1],pos[2] );
996                                                 glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
997                                         } // */
998                                         if((mpControl->mDebugCompavScale>0.) && (TESTGET_FORCE(lev,i,j,k).compAvWeight>0.)) {
999                                                 const float scale = mpControl->mDebugCompavScale*1.;
1000                                                 float dx = TESTGET_FORCE(lev,i,j,k).compAv[0];
1001                                                 float dy = TESTGET_FORCE(lev,i,j,k).compAv[1];
1002                                                 float dz = TESTGET_FORCE(lev,i,j,k).compAv[2];
1003                                                 dx *= scale; dy *= scale; dz *= scale;
1004                                                 glColor3f( 0.2,0.2,1 );
1005                                                 glVertex3f( pos[0],pos[1],pos[2] );
1006                                                 glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
1007                                         } // */
1008                         } // att,maxd
1009                         }
1010                         glEnd();
1011                 }
1012         } // cpssi
1013
1014         //fprintf(stderr,"BLA\n");
1015         glEnable( GL_LIGHTING ); // dont light lines
1016         glShadeModel(GL_SMOOTH);
1017 }
1018
1019 #else // LBM_USE_GUI==1
1020 void LbmFsgrSolver::cpDebugDisplay(int dispset) { }
1021 #endif // LBM_USE_GUI==1
1022
1023