intern/elbeem/intern/solver_control.cpp

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