Blender V2.61 - r43446

solver_control.cpp

Go to the documentation of this file.
00001 
00004 /******************************************************************************
00005  *
00006  * El'Beem - the visual lattice boltzmann freesurface simulator
00007  * All code distributed as part of El'Beem is covered by the version 2 of the 
00008  * GNU General Public License. See the file COPYING for details.  
00009  *
00010  * Copyright 2003-2008 Nils Thuerey 
00011  *
00012  * control extensions
00013  *
00014  *****************************************************************************/
00015 #include "solver_class.h"
00016 #include "solver_relax.h"
00017 #include "particletracer.h"
00018 
00019 #include "solver_control.h"
00020 
00021 #include "controlparticles.h"
00022 
00023 #include "elbeem.h"
00024 
00025 #include "ntl_geometrymodel.h"
00026 
00027 /******************************************************************************
00028  * LbmControlData control set
00029  *****************************************************************************/
00030 
00031 LbmControlSet::LbmControlSet() :
00032     mCparts(NULL), mCpmotion(NULL), mContrPartFile(""), mCpmotionFile(""),
00033     mcForceAtt(0.), mcForceVel(0.), mcForceMaxd(0.), 
00034     mcRadiusAtt(0.), mcRadiusVel(0.), mcRadiusMind(0.), mcRadiusMaxd(0.), 
00035     mcCpScale(1.), mcCpOffset(0.)
00036 {
00037 }
00038 LbmControlSet::~LbmControlSet() {
00039     if(mCparts) delete mCparts;
00040     if(mCpmotion) delete mCpmotion;
00041 }
00042 void LbmControlSet::initCparts() {
00043     mCparts = new ControlParticles();
00044     mCpmotion = new ControlParticles();
00045 }
00046 
00047 
00048 
00049 /******************************************************************************
00050  * LbmControlData control
00051  *****************************************************************************/
00052 
00053 LbmControlData::LbmControlData() : 
00054     mSetForceStrength(0.),
00055     mCons(), 
00056     mCpUpdateInterval(8), // DG: was 16 --> causes problems (big sphere after some time), unstable
00057     mCpOutfile(""), 
00058     mCpForces(), mCpKernel(), mMdKernel(),
00059     mDiffVelCon(1.),
00060     mDebugCpscale(0.), 
00061     mDebugVelScale(0.), 
00062     mDebugCompavScale(0.), 
00063     mDebugAttScale(0.), 
00064     mDebugMaxdScale(0.),
00065     mDebugAvgVelScale(0.)
00066 {
00067 }
00068 
00069 LbmControlData::~LbmControlData() 
00070 {
00071     while (!mCons.empty()) {
00072         delete mCons.back();  mCons.pop_back();
00073     }
00074 }
00075 
00076 
00077 void LbmControlData::parseControldataAttrList(AttributeList *attr) {
00078     // controlpart vars
00079     mSetForceStrength = attr->readFloat("tforcestrength", mSetForceStrength,"LbmControlData", "mSetForceStrength", false);
00080     //errMsg("tforcestrength set to "," "<<mSetForceStrength);
00081     mCpUpdateInterval = attr->readInt("controlparticle_updateinterval", mCpUpdateInterval,"LbmControlData","mCpUpdateInterval", false);
00082     // tracer output file
00083     mCpOutfile = attr->readString("controlparticle_outfile",mCpOutfile,"LbmControlData","mCpOutfile", false);
00084     if(getenv("ELBEEM_CPOUTFILE")) {
00085         string outfile(getenv("ELBEEM_CPOUTFILE"));
00086         mCpOutfile = outfile;
00087         debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPOUTFILE to set mCpOutfile to "<<outfile<<","<<mCpOutfile,7);
00088     }
00089 
00090     for(int cpii=0; cpii<10; cpii++) {
00091         string suffix("");
00092         //if(cpii>0)
00093         {  suffix = string("0"); suffix[0]+=cpii; }
00094         LbmControlSet *cset;
00095         cset = new LbmControlSet();
00096         cset->initCparts();
00097 
00098         cset->mContrPartFile = attr->readString("controlparticle"+suffix+"_file",cset->mContrPartFile,"LbmControlData","cset->mContrPartFile", false);
00099         if((cpii==0) && (getenv("ELBEEM_CPINFILE")) ) {
00100             string infile(getenv("ELBEEM_CPINFILE"));
00101             cset->mContrPartFile = infile;
00102             debMsgStd("LbmControlData::parseAttrList",DM_NOTIFY,"Using envvar ELBEEM_CPINFILE to set mContrPartFile to "<<infile<<","<<cset->mContrPartFile,7);
00103         }
00104 
00105         LbmFloat cpvort=0.;
00106         cset->mcRadiusAtt =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusatt", 0., "LbmControlData","mcRadiusAtt" );
00107         cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
00108         cset->mcRadiusVel =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusvel", 0., "LbmControlData","mcRadiusVel" );
00109         cset->mCparts->setRadiusAtt(cset->mcRadiusAtt.get(0.));
00110         cset->mCparts->setRadiusVel(cset->mcRadiusVel.get(0.));
00111 
00112         // WARNING currently only for first set
00113         //if(cpii==0) {
00114         cset->mcForceAtt  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_attraction", 0. , "LbmControlData","cset->mcForceAtt", false);
00115         cset->mcForceVel  =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_velocity",   0. , "LbmControlData","mcForceVel", false);
00116         cset->mcForceMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_maxdist",    0. , "LbmControlData","mcForceMaxd", false);
00117         cset->mCparts->setInfluenceAttraction(cset->mcForceAtt.get(0.) );
00118         // warning - stores temprorarily, value converted to dt dep. factor
00119         cset->mCparts->setInfluenceVelocity(cset->mcForceVel.get(0.) , 0.01 ); // dummy dt
00120         cset->mCparts->setInfluenceMaxdist(cset->mcForceMaxd.get(0.) );
00121         cpvort =  attr->readFloat("controlparticle"+suffix+"_vorticity",   cpvort, "LbmControlData","cpvort", false);
00122         cset->mCparts->setInfluenceTangential(cpvort);
00123             
00124         cset->mcRadiusMind =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmin", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMind", false);
00125         cset->mcRadiusMaxd =  attr->readChannelSinglePrecFloat("controlparticle"+suffix+"_radiusmax", cset->mcRadiusMind.get(0.), "LbmControlData","mcRadiusMaxd", false);
00126         cset->mCparts->setRadiusMinMaxd(cset->mcRadiusMind.get(0.));
00127         cset->mCparts->setRadiusMaxd(cset->mcRadiusMaxd.get(0.));
00128         //}
00129 
00130         // now local...
00131         //LbmVec cpOffset(0.), cpScale(1.);
00132         LbmFloat cpTimescale = 1.;
00133         string cpMirroring("");
00134 
00135         //cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
00136         //cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
00137         cset->mcCpOffset = attr->readChannelVec3f("controlparticle"+suffix+"_offset", ntlVec3f(0.),"LbmControlData","mcCpOffset", false);
00138         cset->mcCpScale =  attr->readChannelVec3f("controlparticle"+suffix+"_scale",  ntlVec3f(1.), "LbmControlData","mcCpScale", false);
00139         cpTimescale =  attr->readFloat("controlparticle"+suffix+"_timescale",  cpTimescale, "LbmControlData","cpTimescale", false);
00140         cpMirroring =  attr->readString("controlparticle"+suffix+"_mirror",  cpMirroring, "LbmControlData","cpMirroring", false);
00141 
00142         LbmFloat cpsWidth = cset->mCparts->getCPSWith();
00143         cpsWidth =  attr->readFloat("controlparticle"+suffix+"_cpswidth",  cpsWidth, "LbmControlData","cpsWidth", false);
00144         LbmFloat cpsDt = cset->mCparts->getCPSTimestep();
00145         cpsDt =  attr->readFloat("controlparticle"+suffix+"_cpstimestep",  cpsDt, "LbmControlData","cpsDt", false);
00146         LbmFloat cpsTstart = cset->mCparts->getCPSTimeStart();
00147         cpsTstart =  attr->readFloat("controlparticle"+suffix+"_cpststart",  cpsTstart, "LbmControlData","cpsTstart", false);
00148         LbmFloat cpsTend = cset->mCparts->getCPSTimeEnd();
00149         cpsTend =  attr->readFloat("controlparticle"+suffix+"_cpstend",  cpsTend, "LbmControlData","cpsTend", false);
00150         LbmFloat cpsMvmfac = cset->mCparts->getCPSMvmWeightFac();
00151         cpsMvmfac =  attr->readFloat("controlparticle"+suffix+"_cpsmvmfac",  cpsMvmfac, "LbmControlData","cpsMvmfac", false);
00152         cset->mCparts->setCPSWith(cpsWidth);
00153         cset->mCparts->setCPSTimestep(cpsDt);
00154         cset->mCparts->setCPSTimeStart(cpsTstart);
00155         cset->mCparts->setCPSTimeEnd(cpsTend);
00156         cset->mCparts->setCPSMvmWeightFac(cpsMvmfac);
00157 
00158         cset->mCparts->setOffset( vec2L(cset->mcCpOffset.get(0.)) );
00159         cset->mCparts->setScale( vec2L(cset->mcCpScale.get(0.)) );
00160         cset->mCparts->setInitTimeScale( cpTimescale );
00161         cset->mCparts->setInitMirror( cpMirroring );
00162 
00163         int mDebugInit = 0;
00164         mDebugInit = attr->readInt("controlparticle"+suffix+"_debuginit", mDebugInit,"LbmControlData","mDebugInit", false);
00165         cset->mCparts->setDebugInit(mDebugInit);
00166 
00167         // motion particle settings
00168         LbmVec mcpOffset(0.), mcpScale(1.);
00169         LbmFloat mcpTimescale = 1.;
00170         string mcpMirroring("");
00171 
00172         cset->mCpmotionFile = attr->readString("cpmotion"+suffix+"_file",cset->mCpmotionFile,"LbmControlData","mCpmotionFile", false);
00173         mcpTimescale =  attr->readFloat("cpmotion"+suffix+"_timescale",  mcpTimescale, "LbmControlData","mcpTimescale", false);
00174         mcpMirroring =  attr->readString("cpmotion"+suffix+"_mirror",  mcpMirroring, "LbmControlData","mcpMirroring", false);
00175         mcpOffset = vec2L( attr->readVec3d("cpmotion"+suffix+"_offset", vec2P(mcpOffset),"LbmControlData","cpOffset", false) );
00176         mcpScale =  vec2L( attr->readVec3d("cpmotion"+suffix+"_scale",  vec2P(mcpScale), "LbmControlData","cpScale", false) );
00177 
00178         cset->mCpmotion->setOffset( vec2L(mcpOffset) );
00179         cset->mCpmotion->setScale( vec2L(mcpScale) );
00180         cset->mCpmotion->setInitTimeScale( mcpTimescale );
00181         cset->mCpmotion->setInitMirror( mcpMirroring );
00182 
00183         if(cset->mContrPartFile.length()>1) {
00184             errMsg("LbmControlData","Using control particle set "<<cpii<<" file:"<<cset->mContrPartFile<<" cpmfile:"<<cset->mCpmotionFile<<" mirr:'"<<cset->mCpmotion->getInitMirror()<<"' " );
00185             mCons.push_back( cset );
00186         } else {
00187             delete cset;
00188         }
00189     }
00190 
00191     // debug, testing - make sure theres at least an empty set
00192     if(mCons.size()<1) {
00193         mCons.push_back( new LbmControlSet() );
00194         mCons[0]->initCparts();
00195     }
00196 
00197     // take from first set
00198     for(int cpii=1; cpii<(int)mCons.size(); cpii++) {
00199         mCons[cpii]->mCparts->setRadiusMinMaxd( mCons[0]->mCparts->getRadiusMinMaxd() );
00200         mCons[cpii]->mCparts->setRadiusMaxd(    mCons[0]->mCparts->getRadiusMaxd() );
00201         mCons[cpii]->mCparts->setInfluenceAttraction( mCons[0]->mCparts->getInfluenceAttraction() );
00202         mCons[cpii]->mCparts->setInfluenceTangential( mCons[0]->mCparts->getInfluenceTangential() );
00203         mCons[cpii]->mCparts->setInfluenceVelocity( mCons[0]->mCparts->getInfluenceVelocity() , 0.01 ); // dummy dt
00204         mCons[cpii]->mCparts->setInfluenceMaxdist( mCons[0]->mCparts->getInfluenceMaxdist() );
00205     }
00206     
00207     // invert for usage in relax macro
00208     mDiffVelCon =  1.-attr->readFloat("cpdiffvelcon",  mDiffVelCon, "LbmControlData","mDiffVelCon", false);
00209 
00210     mDebugCpscale     =  attr->readFloat("cpdebug_cpscale",  mDebugCpscale, "LbmControlData","mDebugCpscale", false);
00211     mDebugMaxdScale   =  attr->readFloat("cpdebug_maxdscale",  mDebugMaxdScale, "LbmControlData","mDebugMaxdScale", false);
00212     mDebugAttScale    =  attr->readFloat("cpdebug_attscale",  mDebugAttScale, "LbmControlData","mDebugAttScale", false);
00213     mDebugVelScale    =  attr->readFloat("cpdebug_velscale",  mDebugVelScale, "LbmControlData","mDebugVelScale", false);
00214     mDebugCompavScale =  attr->readFloat("cpdebug_compavscale",  mDebugCompavScale, "LbmControlData","mDebugCompavScale", false);
00215     mDebugAvgVelScale =  attr->readFloat("cpdebug_avgvelsc",  mDebugAvgVelScale, "LbmControlData","mDebugAvgVelScale", false);
00216 }
00217 
00218 
00219 void 
00220 LbmFsgrSolver::initCpdata()
00221 {
00222     // enable for cps via env. vars
00223     //if( (getenv("ELBEEM_CPINFILE")) || (getenv("ELBEEM_CPOUTFILE")) ){ mUseTestdata=1; }
00224 
00225     
00226     // manually switch on! if this is zero, nothing is done...
00227     mpControl->mSetForceStrength = this->mTForceStrength = 1.;
00228     mpControl->mCons.clear();
00229     
00230     // init all control fluid objects
00231     int numobjs = (int)(mpGiObjects->size());
00232     for(int o=0; o<numobjs; o++) {
00233         ntlGeometryObjModel *obj = (ntlGeometryObjModel *)(*mpGiObjects)[o];
00234         if(obj->getGeoInitType() & FGI_CONTROL) {
00235             // add new control set per object
00236             LbmControlSet *cset;
00237 
00238             cset = new LbmControlSet();
00239             cset->initCparts();
00240     
00241             // dont load any file
00242             cset->mContrPartFile = string("");
00243 
00244             cset->mcForceAtt = obj->getCpsAttrFStr();
00245             cset->mcRadiusAtt = obj->getCpsAttrFRad();
00246             cset->mcForceVel = obj->getCpsVelFStr();
00247             cset->mcRadiusVel = obj->getCpsVelFRad();
00248 
00249             cset->mCparts->setCPSTimeStart(obj->getCpsTimeStart());
00250             cset->mCparts->setCPSTimeEnd(obj->getCpsTimeEnd());
00251             
00252             if(obj->getCpsQuality() > LBM_EPSILON)
00253                 cset->mCparts->setCPSWith(1.0 / obj->getCpsQuality());
00254             
00255             // this value can be left at 0.5:
00256             cset->mCparts->setCPSMvmWeightFac(0.5);
00257 
00258             mpControl->mCons.push_back( cset );
00259             mpControl->mCons[mpControl->mCons.size()-1]->mCparts->initFromObject(obj);
00260         }
00261     }
00262     
00263     // NT blender integration manual test setup
00264     if(0) {
00265         // manually switch on! if this is zero, nothing is done...
00266         mpControl->mSetForceStrength = this->mTForceStrength = 1.;
00267         mpControl->mCons.clear();
00268 
00269         // add new set
00270         LbmControlSet *cset;
00271 
00272         cset = new LbmControlSet();
00273         cset->initCparts();
00274         // dont load any file
00275         cset->mContrPartFile = string("");
00276 
00277         // set radii for attraction & velocity forces
00278         // set strength of the forces
00279         // don't set directly! but use channels: 
00280         // mcForceAtt, mcForceVel, mcForceMaxd, mcRadiusAtt, mcRadiusVel, mcRadiusMind, mcRadiusMaxd etc.
00281 
00282         // wrong: cset->mCparts->setInfluenceAttraction(1.15); cset->mCparts->setRadiusAtt(1.5);
00283         // right, e.g., to init some constant values:
00284         cset->mcForceAtt = AnimChannel<float>(0.2);
00285         cset->mcRadiusAtt = AnimChannel<float>(0.75);
00286         cset->mcForceVel = AnimChannel<float>(0.2);
00287         cset->mcRadiusVel = AnimChannel<float>(0.75);
00288 
00289         // this value can be left at 0.5:
00290         cset->mCparts->setCPSMvmWeightFac(0.5);
00291 
00292         mpControl->mCons.push_back( cset );
00293 
00294         // instead of reading from file (cset->mContrPartFile), manually init some particles
00295         mpControl->mCons[0]->mCparts->initBlenderTest();
00296 
00297         // other values that might be interesting to change:
00298         //cset->mCparts->setCPSTimestep(0.02);
00299         //cset->mCparts->setCPSTimeStart(0.);
00300         //cset->mCparts->setCPSTimeEnd(1.); 
00301 
00302         //mpControl->mDiffVelCon = 1.; // more rigid velocity control, 0 (default) allows more turbulence
00303     }
00304 
00305     // control particle -------------------------------------------------------------------------------------
00306 
00307     // init cppf stage, use set 0!
00308     if(mCppfStage>0) {
00309         if(mpControl->mCpOutfile.length()<1) mpControl->mCpOutfile = string("cpout"); // use getOutFilename !?
00310         char strbuf[100];
00311         const char *cpFormat = "_d%dcppf%d";
00312 
00313         // initial coarse stage, no input
00314         if(mCppfStage==1) {
00315             mpControl->mCons[0]->mContrPartFile = "";
00316         } else {
00317             // read from prev stage
00318             snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage-1);
00319             mpControl->mCons[0]->mContrPartFile = mpControl->mCpOutfile;
00320             mpControl->mCons[0]->mContrPartFile += strbuf;
00321             mpControl->mCons[0]->mContrPartFile += ".cpart2";
00322         }
00323 
00324         snprintf(strbuf,100, cpFormat  ,LBMDIM,mCppfStage);
00325         mpControl->mCpOutfile += strbuf;
00326     } // */
00327     
00328     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00329         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00330         ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
00331 
00332         // now set with real dt
00333         cparts->setInfluenceVelocity( mpControl->mCons[cpssi]->mcForceVel.get(0.), mLevel[mMaxRefine].timestep);
00334         cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
00335         cparts->setCharLength( mLevel[mMaxRefine].nodeSize );
00336         errMsg("LbmControlData","CppfStage "<<mCppfStage<<" in:"<<mpControl->mCons[cpssi]->mContrPartFile<<
00337                 " out:"<<mpControl->mCpOutfile<<" cl:"<< cparts->getCharLength() );
00338 
00339         // control particle test init
00340         if(mpControl->mCons[cpssi]->mCpmotionFile.length()>=1) cpmotion->initFromTextFile(mpControl->mCons[cpssi]->mCpmotionFile);
00341         // not really necessary...
00342         //? cparts->setFluidSpacing( mLevel[mMaxRefine].nodeSize    ); // use grid coords!?
00343         //? cparts->calculateKernelWeight();
00344         //? debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - motion inited: "<<cparts->getSize() ,10);
00345 
00346         // ensure both are on for env. var settings
00347         // when no particles, but outfile enabled, initialize
00348         const int lev = mMaxRefine;
00349         if((mpParticles) && (mpControl->mCpOutfile.length()>=1) && (cpssi==0)) {
00350             // check if auto num
00351             if( (mpParticles->getNumInitialParticles()<=1) && 
00352                     (mpParticles->getNumParticles()<=1) ) { // initParticles done afterwards anyway
00353                 int tracers = 0;
00354                 const int workSet = mLevel[lev].setCurr;
00355                 FSGR_FORIJK_BOUNDS(lev) { 
00356                     if(RFLAG(lev,i,j,k, workSet)&(CFFluid)) tracers++;
00357                 }
00358                 if(LBMDIM==3) tracers /= 8;
00359                 else          tracers /= 4;
00360                 mpParticles->setNumInitialParticles(tracers);
00361                 mpParticles->setDumpTextFile(mpControl->mCpOutfile);
00362                 debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - set tracers #"<<tracers<<", actual #"<<mpParticles->getNumParticles() ,10);
00363             }
00364             if(mpParticles->getDumpTextInterval()<=0.) {
00365                 mpParticles->setDumpTextInterval(mLevel[lev].timestep * mLevel[lev].lSizex);
00366                 debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles - dump delta t not set, using dti="<< mpParticles->getDumpTextInterval()<<", sim dt="<<mLevel[lev].timestep, 5 );
00367             }
00368             mpParticles->setDumpParts(true); // DEBUG? also dump as particle system
00369         }
00370 
00371         if(mpControl->mCons[cpssi]->mContrPartFile.length()>=1) cparts->initFromTextFile(mpControl->mCons[cpssi]->mContrPartFile);
00372         cparts->setFluidSpacing( mLevel[lev].nodeSize   ); // use grid coords!?
00373         cparts->calculateKernelWeight();
00374         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);
00375     } // cpssi
00376 
00377     if(getenv("ELBEEM_CPINFILE")) {
00378         this->mTForceStrength = 1.0;
00379     }
00380     this->mTForceStrength = mpControl->mSetForceStrength;
00381     if(mpControl->mCpOutfile.length()>=1) mpParticles->setDumpTextFile(mpControl->mCpOutfile);
00382 
00383     // control particle  init end -------------------------------------------------------------------------------------
00384 
00385     // make sure equiv to solver init
00386     if(this->mTForceStrength>0.) { \
00387         mpControl->mCpForces.resize( mMaxRefine+1 );
00388         for(int lev = 0; lev<=mMaxRefine; lev++) {
00389             LONGINT rcellSize = (mLevel[lev].lSizex*mLevel[lev].lSizey*mLevel[lev].lSizez);
00390             debMsgStd("LbmFsgrSolver::initControl",DM_MSG,"mCpForces init, lev="<<lev<<" rcs:"<<(int)(rcellSize+4)<<","<<(rcellSize*sizeof(ControlForces)/(1024*1024)), 9 );
00391             mpControl->mCpForces[lev].resize( (int)(rcellSize+4) );
00392             //for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces.push_back( ControlForces() );
00393             for(int i=0 ;i<rcellSize; i++) mpControl->mCpForces[lev][i].resetForces();
00394         }
00395     } // on?
00396 
00397     debMsgStd("LbmFsgrSolver::initCpdata",DM_MSG,"ControlParticles #mCons "<<mpControl->mCons.size()<<" done", 6);
00398 }
00399 
00400 
00401 #define CPODEBUG 0
00402 //define CPINTER ((int)(mpControl->mCpUpdateInterval))
00403 
00404 #define KERN(x,y,z)    mpControl->mCpKernel[ (((z)*cpkarWidth + (y))*cpkarWidth + (x)) ]
00405 #define MDKERN(x,y,z)  mpControl->mMdKernel[ (((z)*mdkarWidth + (y))*mdkarWidth + (x)) ]
00406 
00407 #define BOUNDCHECK(x,low,high)  ( ((x)<low) ? low : (((x)>high) ? high : (x) )  )
00408 #define BOUNDSKIP(x,low,high)  ( ((x)<low) || ((x)>high) )
00409 
00410 void 
00411 LbmFsgrSolver::handleCpdata()
00412 {
00413     myTime_t cpstart = getTime();
00414     int cpChecks=0;
00415     int cpInfs=0;
00416     //debMsgStd("ControlData::handleCpdata",DM_MSG,"called... "<<this->mTForceStrength,1);
00417 
00418     // add cp influence
00419   if((true) && (this->mTForceStrength>0.)) {
00420         // ok continue...
00421     } // on off
00422     else {
00423         return;
00424     }
00425     
00426     // check if we have control objects
00427     if(mpControl->mCons.size()==0)
00428         return;
00429     
00430     if((mpControl->mCpUpdateInterval<1) || (this->mStepCnt%mpControl->mCpUpdateInterval==0)) {
00431         // do full reinit later on... 
00432     }
00433     else if(this->mStepCnt>mpControl->mCpUpdateInterval) {
00434         // only reinit new cells
00435         // TODO !? remove loop dependance!?
00436 #define NOFORCEENTRY(lev, i,j,k) (LBMGET_FORCE(lev, i,j,k).maxDistance==CPF_MAXDINIT) 
00437         // interpolate missing
00438         for(int lev=0; lev<=mMaxRefine; lev++) {
00439         FSGR_FORIJK_BOUNDS(lev) { 
00440             if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFFluid|CFInter) )  
00441             //if( (RFLAG(lev,i,j,k, mLevel[lev].setCurr)) & (CFInter) )  
00442             //if(0)
00443             { // only check new inter? RFLAG?check
00444                 if(NOFORCEENTRY(lev, i,j,k)) {
00445                     //errMsg("CP","FE_MISSING at "<<PRINT_IJK<<" f"<<LBMGET_FORCE(lev, i,j,k).weightAtt<<" md"<<LBMGET_FORCE(lev, i,j,k).maxDistance );
00446 
00447                     LbmFloat nbs=0.;
00448                     ControlForces vals;
00449                     vals.resetForces(); vals.maxDistance = 0.;
00450                     for(int l=1; l<this->cDirNum; l++) { 
00451                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
00452                         //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 );
00453                         if(!NOFORCEENTRY(lev, ni,nj,nk)) {
00454                             //? vals.weightAtt   += LBMGET_FORCE(lev, ni,nj,nk).weightAtt;
00455                             //? vals.forceAtt    += LBMGET_FORCE(lev, ni,nj,nk).forceAtt;
00456                             vals.maxDistance += LBMGET_FORCE(lev, ni,nj,nk).maxDistance;
00457                             vals.forceMaxd   += LBMGET_FORCE(lev, ni,nj,nk).forceMaxd;
00458                             vals.weightVel   += LBMGET_FORCE(lev, ni,nj,nk).weightVel;
00459                             vals.forceVel    += LBMGET_FORCE(lev, ni,nj,nk).forceVel; 
00460                             // ignore att/compAv/avgVel here for now
00461                             nbs += 1.;
00462                         }
00463                     }
00464                     if(nbs>0.) {
00465                         nbs = 1./nbs;
00466                         //? LBMGET_FORCE(lev, i,j,k).weightAtt   =  vals.weightAtt*nbs;
00467                         //? LBMGET_FORCE(lev, i,j,k).forceAtt    = vals.forceAtt*nbs;
00468                         LBMGET_FORCE(lev, i,j,k).maxDistance = vals.maxDistance*nbs;
00469                         LBMGET_FORCE(lev, i,j,k).forceMaxd   =  vals.forceMaxd*nbs;
00470                         LBMGET_FORCE(lev, i,j,k).weightVel   =  vals.weightVel*nbs;
00471                         LBMGET_FORCE(lev, i,j,k).forceVel    = vals.forceVel*nbs;
00472                     }
00473                         /*ControlForces *ff = &LBMGET_FORCE(lev, i,j,k);  // DEBUG
00474                         errMsg("CP","FE_MISSING rec at "<<PRINT_IJK // DEBUG
00475                                 <<" w:"<<ff->weightAtt<<" wa:" <<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
00476                                 <<" v:"<<ff->weightVel<<" wv:" <<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
00477                                 <<" v:"<<ff->maxDistance<<" wv:" <<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  ); // DEBUG */
00478                     // else errMsg("CP","FE_MISSING rec at "<<PRINT_IJK<<" failed!"); // DEBUG
00479                     
00480                 }
00481             }
00482         }} // ijk, lev
00483 
00484         // mStepCnt > mpControl->mCpUpdateInterval
00485         return;
00486     } else {
00487         // nothing to do ...
00488         return;
00489     }
00490 
00491     // reset
00492     for(int lev=0; lev<=mMaxRefine; lev++) {
00493         FSGR_FORIJK_BOUNDS(lev) { LBMGET_FORCE(lev,i,j,k).resetForces(); }
00494     }
00495     // do setup for coarsest level
00496     const int coarseLev = 0;
00497     const int fineLev = mMaxRefine;
00498 
00499     // init for current time
00500     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00501         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00502         LbmControlSet *cset = mpControl->mCons[cpssi];
00503         
00504         cparts->setRadiusAtt(cset->mcRadiusAtt.get(mSimulationTime));
00505         cparts->setRadiusVel(cset->mcRadiusVel.get(mSimulationTime));
00506         cparts->setInfluenceAttraction(cset->mcForceAtt.get(mSimulationTime) );
00507         cparts->setInfluenceMaxdist(cset->mcForceMaxd.get(mSimulationTime) );
00508         cparts->setRadiusMinMaxd(cset->mcRadiusMind.get(mSimulationTime));
00509         cparts->setRadiusMaxd(cset->mcRadiusMaxd.get(mSimulationTime));
00510         cparts->calculateKernelWeight(); // always necessary!?
00511         cparts->setOffset( vec2L(cset->mcCpOffset.get(mSimulationTime)) );
00512         cparts->setScale(  vec2L(cset->mcCpScale.get(mSimulationTime)) );
00513 
00514         cparts->setInfluenceVelocity( cset->mcForceVel.get(mSimulationTime), mLevel[fineLev].timestep );
00515         cparts->setLastOffset( vec2L(cset->mcCpOffset.get(mSimulationTime-mLevel[fineLev].timestep)) );
00516         cparts->setLastScale(  vec2L(cset->mcCpScale.get(mSimulationTime-mLevel[fineLev].timestep)) );
00517         
00518     }
00519 
00520     // check actual values
00521     LbmFloat iatt  = ABS(mpControl->mCons[0]->mCparts->getInfluenceAttraction());
00522     LbmFloat ivel  = mpControl->mCons[0]->mCparts->getInfluenceVelocity();
00523     LbmFloat imaxd = mpControl->mCons[0]->mCparts->getInfluenceMaxdist();
00524     //errMsg("FINCIT","iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
00525     for(int cpssi=1; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00526         LbmFloat iatt2  = ABS(mpControl->mCons[cpssi]->mCparts->getInfluenceAttraction());
00527         LbmFloat ivel2  = mpControl->mCons[cpssi]->mCparts->getInfluenceVelocity();
00528         LbmFloat imaxd2 = mpControl->mCons[cpssi]->mCparts->getInfluenceMaxdist();
00529         
00530         // we allow negative attraction force here!
00531         if(iatt2 > iatt)   iatt = iatt2;
00532         
00533         if(ivel2 >ivel)   ivel = ivel2;
00534         if(imaxd2>imaxd)  imaxd= imaxd2;
00535         //errMsg("FINCIT"," "<<cpssi<<" iatt2="<<iatt2<<" ivel2="<<ivel2<<" imaxd2="<<imaxd<<" NEW "<<" iatt="<<iatt<<" ivel="<<ivel<<" imaxd="<<imaxd);
00536     }
00537 
00538     if(iatt==0. && ivel==0. && imaxd==0.) {
00539         debMsgStd("ControlData::initControl",DM_MSG,"Skipped, all zero...",4);
00540         return;
00541     }
00542     //iatt  = mpControl->mCons[1]->mCparts->getInfluenceAttraction(); //ivel  = mpControl->mCons[1]->mCparts->getInfluenceVelocity(); //imaxd = mpControl->mCons[1]->mCparts->getInfluenceMaxdist(); // TTTTTT
00543 
00544     // do control setup
00545     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00546         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00547         ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
00548 
00549         // TEST!?
00550         bool radmod = false;
00551         const LbmFloat minRadSize = mLevel[coarseLev].nodeSize * 1.5;
00552         if((cparts->getRadiusAtt()>0.) && (cparts->getRadiusAtt()<minRadSize) && (!radmod) ) {
00553             LbmFloat radfac = minRadSize / cparts->getRadiusAtt(); radmod=true;
00554             debMsgStd("ControlData::initControl",DM_MSG,"Modified radii att, fac="<<radfac, 7);
00555             cparts->setRadiusAtt(cparts->getRadiusAtt()*radfac);
00556             cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
00557             cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
00558             cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
00559         } else if((cparts->getRadiusVel()>0.) && (cparts->getRadiusVel()<minRadSize) && (!radmod) ) {
00560             LbmFloat radfac = minRadSize / cparts->getRadiusVel();
00561             debMsgStd("ControlData::initControl",DM_MSG,"Modified radii vel, fac="<<radfac, 7);
00562             cparts->setRadiusVel(cparts->getRadiusVel()*radfac);
00563             cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
00564             cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
00565         } else if((cparts->getRadiusMaxd()>0.) && (cparts->getRadiusMaxd()<minRadSize) && (!radmod) ) {
00566             LbmFloat radfac = minRadSize / cparts->getRadiusMaxd();
00567             debMsgStd("ControlData::initControl",DM_MSG,"Modified radii maxd, fac="<<radfac, 7);
00568             cparts->setRadiusMaxd(cparts->getRadiusMaxd()*radfac);
00569             cparts->setRadiusMinMaxd(cparts->getRadiusMinMaxd()*radfac);
00570         }
00571         if(radmod) {
00572             debMsgStd("ControlData::initControl",DM_MSG,"Modified radii: att="<< 
00573                 cparts->getRadiusAtt()<<", vel=" << cparts->getRadiusVel()<<", maxd=" <<
00574                 cparts->getRadiusMaxd()<<", mind=" << cparts->getRadiusMinMaxd() ,5);
00575         }
00576 
00577         cpmotion->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), NULL );
00578         cparts->prepareControl( mSimulationTime+((LbmFloat)mpControl->mCpUpdateInterval)*(mpParam->getTimestep()), mpParam->getTimestep(), cpmotion );
00579     }
00580 
00581     // do control...
00582     for(int lev=0; lev<=mMaxRefine; lev++) {
00583         LbmFloat levVolume = 1.;
00584         LbmFloat levForceScale = 1.;
00585         for(int ll=lev; ll<mMaxRefine; ll++) {
00586             if(LBMDIM==3) levVolume *= 8.;
00587             else levVolume *= 4.;
00588             levForceScale *= 2.;
00589         }
00590         errMsg("LbmFsgrSolver::handleCpdata","levVolume="<<levVolume<<" levForceScale="<<levForceScale );
00591     //todo: scale velocity, att by level timestep!?
00592         
00593     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00594         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00595         // ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
00596         
00597         // if control set is not active skip it
00598         if((cparts->getControlTimStart() > mSimulationTime) || (cparts->getControlTimEnd() < mLastSimTime))
00599         {
00600             continue;
00601         }
00602 
00603         const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
00604         LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
00605         LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
00606         LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
00607 #if LBMDIM==2
00608         gsz = gsx;
00609 #endif
00610         LbmFloat goffx = mvGeoStart[0];
00611         LbmFloat goffy = mvGeoStart[1];
00612         LbmFloat goffz = mvGeoStart[2];
00613 
00614         //const LbmFloat cpwIncFac = 2.0;
00615         // max to two thirds of domain size
00616         const int cpw = MIN( mLevel[lev].lSizex/3, MAX( (int)( cparts->getRadiusAtt()  /gsx) +1  , 2) ); // normal kernel, att,vel
00617         const int cpkarWidth = 2*cpw+1;
00618         mpControl->mCpKernel.resize(cpkarWidth* cpkarWidth* cpkarWidth);
00619         ControlParticle cpt; cpt.reset();
00620         cpt.pos = LbmVec( (gsx*(LbmFloat)cpw)+goffx, (gsy*(LbmFloat)cpw)+goffy, (gsz*(LbmFloat)cpw)+goffz );  // optimize?
00621         cpt.density = 0.5; cpt.densityWeight = 0.5;
00622 #if LBMDIM==3
00623         for(int k= 0; k<cpkarWidth; ++k) {
00624 #else // LBMDIM==3
00625         { int k = cpw;
00626 #endif 
00627             for(int j= 0; j<cpkarWidth; ++j) 
00628                 for(int i= 0; i<cpkarWidth; ++i) {
00629                     KERN(i,j,k).resetForces();
00630                     //LbmFloat dx = i-cpw; LbmFloat dy = j-cpw; LbmFloat dz = k-cpw;
00631                     //LbmVec dv = ( LbmVec(dx,dy,dz) );
00632                     //LbmFloat dl = norm( dv ); //LbmVec dir = dv / dl;
00633                     LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
00634                     cparts->calculateCpInfluenceOpt( &cpt, pos, LbmVec(0,0,0), &KERN(i,j,k)  ,1. );
00635                     /*if((CPODEBUG)&&(k==cpw)) errMsg("kern"," at "<<PRINT_IJK<<" pos"<<pos<<" cpp"<<cpt.pos 
00636                         <<" 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] )  
00637                         <<" 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] )  
00638                         <<" 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] )  ); // */
00639                     KERN(i,j,k).weightAtt *= 2.;
00640                     KERN(i,j,k).forceAtt *= 2.;
00641                     //KERN(i,j,k).forceAtt[1] *= 2.; KERN(i,j,k).forceAtt[2] *= 2.;
00642                     KERN(i,j,k).weightVel *= 2.;
00643                     KERN(i,j,k).forceVel *= 2.;
00644                     //KERN(i,j,k).forceVel[1] *= 2.; KERN(i,j,k).forceVel[2] *= 2.;
00645                 }
00646         }
00647 
00648         if(CPODEBUG) errMsg("cpw"," = "<<cpw<<" f"<< cparts->getRadiusAtt()<<" gsx"<<gsx<<" kpw"<<cpkarWidth); // DEBUG
00649         // first cp loop - add att and vel forces
00650         for(int cppi=0; cppi<cparts->getSize(); cppi++) {
00651             ControlParticle *cp = cparts->getParticle(cppi);
00652             if(cp->influence<=0.) continue;
00653             const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
00654             const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
00655             int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
00656             /*if( ((LBMDIM==3)&&(BOUNDSKIP(cpk - cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
00657                 ((LBMDIM==3)&&(BOUNDSKIP(cpk + cpwsm, getForZMinBnd(), getForZMaxBnd(lev) ))) ||
00658                 BOUNDSKIP(cpj - cpwsm, 0, mLevel[lev].lSizey ) ||
00659                 BOUNDSKIP(cpj + cpwsm, 0, mLevel[lev].lSizey ) ||
00660                 BOUNDSKIP(cpi - cpwsm, 0, mLevel[lev].lSizex ) ||
00661                 BOUNDSKIP(cpi + cpwsm, 0, mLevel[lev].lSizex ) ) {
00662                 continue;
00663                 } // */
00664             int is,ie,js,je,ks,ke;
00665             ks = BOUNDCHECK(cpk - cpw, getForZMinBnd(), getForZMaxBnd(lev) );
00666             ke = BOUNDCHECK(cpk + cpw, getForZMinBnd(), getForZMaxBnd(lev) );
00667             js = BOUNDCHECK(cpj - cpw, 0, mLevel[lev].lSizey );
00668             je = BOUNDCHECK(cpj + cpw, 0, mLevel[lev].lSizey );
00669             is = BOUNDCHECK(cpi - cpw, 0, mLevel[lev].lSizex );
00670             ie = BOUNDCHECK(cpi + cpw, 0, mLevel[lev].lSizex );
00671             if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
00672             if(CPODEBUG) errMsg("cppft","i"<<cppi<<" cpw"<<cpw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
00673             cpInfs++;
00674 
00675             for(int k= ks; k<ke; ++k) {
00676                 for(int j= js; j<je; ++j) {
00677 
00678                     CellFlagType *pflag = &RFLAG(lev,is,j,k, mLevel[lev].setCurr);
00679                     ControlForces *kk = &KERN( is-cpi+cpw, j-cpj+cpw, k-cpk+cpw);
00680                     ControlForces *ff = &LBMGET_FORCE(lev,is,j,k); 
00681                     pflag--; kk--; ff--;
00682 
00683                     for(int i= is; i<ie; ++i) {
00684                         // first cp loop (att,vel)
00685                         pflag++; kk++; ff++;
00686 
00687                         //add weight for bnd cells
00688                         const LbmFloat pwforce = kk->weightAtt;
00689                         // control particle mod,
00690                         // dont add multiple CFFluid fsgr boundaries
00691                         if(lev==mMaxRefine) {
00692                             //if( ( ((*pflag)&(CFFluid )) && (lev==mMaxRefine) ) ||
00693                                     //( ((*pflag)&(CFGrNorm)) && (lev <mMaxRefine) ) ) {
00694                             if((*pflag)&(CFFluid|CFUnused)) {
00695                                 // check not fromcoarse?
00696                                 cp->density += levVolume* kk->weightAtt; // old CFFluid
00697                             } else if( (*pflag) & (CFEmpty) ) {  
00698                                 cp->density -= levVolume* 0.5; 
00699                             } else { //if( ((*pflag) & (CFBnd)) ) {  
00700                                 cp->density -= levVolume* 0.2;  // penalty
00701                             }
00702                         } else {
00703                             //if((*pflag)&(CFGrNorm)) {
00704                                 //cp->density += levVolume* kk->weightAtt; // old CFFluid
00705                             //} 
00706                         }
00707                         //else if(!((*pflag) & (CFUnused)) ) {  cp->density -= levVolume* 0.2; } // penalty
00708 
00709                         if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
00710                         {
00711 
00712                             cpChecks++;
00713                             //const LbmFloat pwforce = kk->weightAtt;
00714                             LbmFloat pwvel = kk->weightVel;
00715                             if((pwforce==0.)&&(pwvel==0.)) { continue; }
00716                             ff->weightAtt += 1e-6; // for distance
00717 
00718                             if(pwforce>0.) {
00719                                 ff->weightAtt += pwforce *cp->densityWeight *cp->influence;
00720                                 ff->forceAtt += kk->forceAtt *levForceScale *cp->densityWeight *cp->influence;
00721 
00722                                 // old fill handling here
00723                                 const int workSet =mLevel[lev].setCurr;
00724                                 LbmFloat ux=0., uy=0., uz=0.;
00725                                 FORDF1{  
00726                                     const LbmFloat dfn = QCELL(lev, i,j,k, workSet, l);
00727                                     ux  += (this->dfDvecX[l]*dfn);
00728                                     uy  += (this->dfDvecY[l]*dfn); 
00729                                     uz  += (this->dfDvecZ[l]*dfn); 
00730                                 }
00731                                 // control particle mod
00732                                 cp->avgVelWeight += levVolume*pwforce;
00733                                 cp->avgVelAcc += LbmVec(ux,uy,uz) * levVolume*pwforce;
00734                             }
00735 
00736                             if(pwvel>0.) {
00737                                 // TODO make switch? vel.influence depends on density weight... 
00738                                 // (reduced lowering with 0.75 factor)
00739                                 pwvel *=  cp->influence *(1.-0.75*cp->densityWeight);
00740                                 // control particle mod
00741                                 // todo use Omega instead!?
00742                                 ff->forceVel += cp->vel*levVolume*pwvel * velLatticeScale; // levVolume?
00743                                 ff->weightVel += levVolume*pwvel; // levVolume?
00744                                 ff->compAv += cp->avgVel*levVolume*pwvel; // levVolume?
00745                                 ff->compAvWeight += levVolume*pwvel; // levVolume?
00746                             }
00747 
00748                             if(CPODEBUG) errMsg("cppft","i"<<cppi<<" at "<<PRINT_IJK<<" kern:"<<
00749                                     PRINT_VEC(i-cpi+cpw, j-cpj+cpw, k-cpk+cpw )
00750                                     //<<" w:"<<ff->weightAtt<<" wa:"
00751                                     //<<PRINT_VEC( ff->forceAtt[0],ff->forceAtt[1],ff->forceAtt[2] )  
00752                                     //<<" v:"<<ff->weightVel<<" wv:"
00753                                     //<<PRINT_VEC( ff->forceVel[0],ff->forceVel[1],ff->forceVel[2] )  
00754                                     //<<" v:"<<ff->maxDistance<<" wv:"
00755                                     //<<PRINT_VEC( ff->forceMaxd[0],ff->forceMaxd[1],ff->forceMaxd[2] )  
00756                                     );
00757                         } // celltype
00758 
00759                     } // ijk
00760                 } // ijk
00761             } // ijk
00762         } // cpi, end first cp loop (att,vel)
00763         debMsgStd("LbmFsgrSolver::handleCpdata",DM_MSG,"Force cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
00764     } //cpssi
00765     } // lev
00766 
00767     // second loop
00768     for(int lev=0; lev<=mMaxRefine; lev++) {
00769         LbmFloat levVolume = 1.;
00770         LbmFloat levForceScale = 1.;
00771         for(int ll=lev; ll<mMaxRefine; ll++) {
00772             if(LBMDIM==3) levVolume *= 8.;
00773             else levVolume *= 4.;
00774             levForceScale *= 2.;
00775         }
00776     // prepare maxd forces
00777     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00778         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00779 
00780         // WARNING copied from above!
00781         const LbmFloat velLatticeScale = mLevel[lev].timestep/mLevel[lev].nodeSize;
00782         LbmFloat gsx = ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex);
00783         LbmFloat gsy = ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey);
00784         LbmFloat gsz = ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez);
00785 #if LBMDIM==2
00786         gsz = gsx;
00787 #endif
00788         LbmFloat goffx = mvGeoStart[0];
00789         LbmFloat goffy = mvGeoStart[1];
00790         LbmFloat goffz = mvGeoStart[2];
00791 
00792         //const LbmFloat cpwIncFac = 2.0;
00793         const int mdw = MIN( mLevel[lev].lSizex/2, MAX( (int)( cparts->getRadiusMaxd() /gsx) +1  , 2) ); // wide kernel, md
00794         const int mdkarWidth = 2*mdw+1;
00795         mpControl->mMdKernel.resize(mdkarWidth* mdkarWidth* mdkarWidth);
00796         ControlParticle cpt; cpt.reset();
00797         cpt.density = 0.5; cpt.densityWeight = 0.5;
00798         cpt.pos = LbmVec( (gsx*(LbmFloat)mdw)+goffx, (gsy*(LbmFloat)mdw)+goffy, (gsz*(LbmFloat)mdw)+goffz );  // optimize?
00799 #if LBMDIM==3
00800         for(int k= 0; k<mdkarWidth; ++k) {
00801 #else // LBMDIM==3
00802         { int k = mdw;
00803 #endif 
00804             for(int j= 0; j<mdkarWidth; ++j) 
00805                 for(int i= 0; i<mdkarWidth; ++i) {
00806                     MDKERN(i,j,k).resetForces();
00807                     LbmVec pos = LbmVec( (gsx*(LbmFloat)i)+goffx, (gsy*(LbmFloat)j)+goffy, (gsz*(LbmFloat)k)+goffz );  // optimize?
00808                     cparts->calculateMaxdForce( &cpt, pos, &MDKERN(i,j,k)  );
00809                 }
00810         }
00811 
00812         // second cpi loop, maxd forces
00813         if(cparts->getInfluenceMaxdist()>0.) {
00814             for(int cppi=0; cppi<cparts->getSize(); cppi++) {
00815                 ControlParticle *cp = cparts->getParticle(cppi);
00816                 if(cp->influence<=0.) continue;
00817                 const int cpi = (int)( (cp->pos[0]-goffx)/gsx );
00818                 const int cpj = (int)( (cp->pos[1]-goffy)/gsy );
00819                 int       cpk = (int)( (cp->pos[2]-goffz)/gsz );
00820 
00821                 int is,ie,js,je,ks,ke;
00822                 ks = BOUNDCHECK(cpk - mdw, getForZMinBnd(), getForZMaxBnd(lev) );
00823                 ke = BOUNDCHECK(cpk + mdw, getForZMinBnd(), getForZMaxBnd(lev) );
00824                 js = BOUNDCHECK(cpj - mdw, 0, mLevel[lev].lSizey );
00825                 je = BOUNDCHECK(cpj + mdw, 0, mLevel[lev].lSizey );
00826                 is = BOUNDCHECK(cpi - mdw, 0, mLevel[lev].lSizex );
00827                 ie = BOUNDCHECK(cpi + mdw, 0, mLevel[lev].lSizex );
00828                 if(LBMDIM==2) { cpk = 0; ks = 0; ke = 1; }
00829                 if(CPODEBUG) errMsg("cppft","i"<<cppi<<" mdw"<<mdw<<" gpos"<<PRINT_VEC(cpi,cpj,cpk)<<"   i:"<<is<<","<<ie<<"   j:"<<js<<","<<je<<"   k:"<<ks<<","<<ke<<" "); // DEBUG
00830                 cpInfs++;
00831 
00832                 for(int k= ks; k<ke; ++k)
00833                  for(int j= js; j<je; ++j) {
00834                     CellFlagType *pflag = &RFLAG(lev,is-1,j,k, mLevel[lev].setCurr);
00835                     for(int i= is; i<ie; ++i) {
00836                         // second cpi loop, maxd forces
00837                         pflag++;
00838                         if( (*pflag) & (CFFluid|CFInter) )  // RFLAG_check
00839                         {
00840                             cpChecks++;
00841                             ControlForces *ff = &LBMGET_FORCE(lev,i,j,k); 
00842                             if(ff->weightAtt == 0.) {
00843                                 ControlForces *kk = &MDKERN( i-cpi+mdw, j-cpj+mdw, k-cpk+mdw);
00844                                 const LbmFloat pmdf = kk->maxDistance;
00845                                 if((ff->maxDistance > pmdf) || (ff->maxDistance<0.))
00846                                     ff->maxDistance = pmdf;
00847                                 ff->forceMaxd = kk->forceMaxd;
00848                                 // todo use Omega instead!?
00849                                 ff->forceVel = cp->vel* velLatticeScale;
00850                             }
00851                         } // celltype
00852                 } } // ijk
00853             } // cpi, md loop 
00854         } // maxd inf>0 */
00855 
00856 
00857         debMsgStd("ControlData::initControl",DM_MSG,"Maxd cpgrid "<<cpssi<<" generated checks:"<<cpChecks<<" infs:"<<cpInfs ,9);
00858     } //cpssi
00859 
00860     // normalize, only done once for the whole array
00861     mpControl->mCons[0]->mCparts->finishControl( mpControl->mCpForces[lev], iatt,ivel,imaxd );
00862 
00863     } // lev loop
00864 
00865     myTime_t cpend = getTime(); 
00866     debMsgStd("ControlData::handleCpdata",DM_MSG,"Time for cpgrid generation:"<< getTimeString(cpend-cpstart)<<", checks:"<<cpChecks<<" infs:"<<cpInfs<<" " ,8);
00867 
00868     // warning, may return  before
00869 }
00870 
00871 #if LBM_USE_GUI==1
00872 
00873 #define USE_GLUTILITIES
00874 #include "../gui/gui_utilities.h"
00875 
00876 void LbmFsgrSolver::cpDebugDisplay(int dispset)
00877 {
00878     for(int cpssi=0; cpssi<(int)mpControl->mCons.size(); cpssi++) {
00879         ControlParticles *cparts = mpControl->mCons[cpssi]->mCparts;
00880         //ControlParticles *cpmotion = mpControl->mCons[cpssi]->mCpmotion;
00881         // display cp parts
00882         const bool cpCubes = false;
00883         const bool cpDots = true;
00884         const bool cpCpdist = true;
00885         const bool cpHideIna = true;
00886         glShadeModel(GL_FLAT);
00887         glDisable( GL_LIGHTING ); // dont light lines
00888 
00889         // dot influence
00890         if((mpControl->mDebugCpscale>0.) && cpDots) {
00891             glPointSize(mpControl->mDebugCpscale * 8.);
00892             glBegin(GL_POINTS);
00893             for(int i=0; i<cparts->getSize(); i++) {
00894                 if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
00895                 ntlVec3Gfx org( vec2G(cparts->getParticle(i)->pos ) );
00896                 //LbmFloat halfsize = 0.5; 
00897                 LbmFloat scale = cparts->getParticle(i)->densityWeight;
00898                 //glColor4f( scale,scale,scale,scale );
00899                 glColor4f( 0.,scale,0.,scale );
00900                 glVertex3f( org[0],org[1],org[2] );
00901                 //errMsg("lbmDebugDisplay","CP "<<i<<" at "<<org); // DEBUG
00902             }
00903             glEnd();
00904         }
00905 
00906         // cp positions
00907         if((mpControl->mDebugCpscale>0.) && cpDots) {
00908             glPointSize(mpControl->mDebugCpscale * 3.);
00909             glBegin(GL_POINTS); 
00910             glColor3f( 0,1,0 );
00911         }
00912         for(int i=0; i<cparts->getSize(); i++) {
00913             if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
00914             ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
00915             LbmFloat halfsize = 0.5; 
00916             LbmFloat scale = cparts->getRadiusAtt() * cparts->getParticle(i)->densityWeight;
00917             if(cpCubes){    glLineWidth( 1 ); 
00918                 glColor3f( 1,1,1 );
00919                 ntlVec3Gfx s = org-(halfsize * (scale)); 
00920                 ntlVec3Gfx e = org+(halfsize * (scale)); 
00921                 drawCubeWire( s,e ); }
00922             if((mpControl->mDebugCpscale>0.) && cpDots) {
00923                 glVertex3f( org[0],org[1],org[2] );
00924             }
00925         }
00926         if(cpDots) glEnd();
00927 
00928         if(mpControl->mDebugAvgVelScale>0.) {
00929             const float scale = mpControl->mDebugAvgVelScale;
00930 
00931             glColor3f( 1.0,1.0,1 );
00932             glBegin(GL_LINES); 
00933             for(int i=0; i<cparts->getSize(); i++) {
00934                 if((cpHideIna)&&( (cparts->getParticle(i)->influence<=0.) || (cparts->getParticle(i)->size<=0.) )) continue;
00935                 ntlVec3Gfx  org( vec2G(cparts->getParticle(i)->pos ) );
00936 
00937                 //errMsg("CPAVGVEL","i"<<i<<" pos"<<org<<" av"<<cparts->getParticle(i)->avgVel);// DEBUG
00938                 float dx = cparts->getParticle(i)->avgVel[0];
00939                 float dy = cparts->getParticle(i)->avgVel[1];
00940                 float dz = cparts->getParticle(i)->avgVel[2];
00941                 dx *= scale; dy *= scale; dz *= scale;
00942                 glVertex3f( org[0],org[1],org[2] );
00943                 glVertex3f( org[0]+dx,org[1]+dy,org[2]+dz );
00944             }
00945             glEnd();
00946         } // */
00947 
00948         if( (LBMDIM==2) && (cpCpdist) ) {
00949             
00950             // debug, for use of e.g. LBMGET_FORCE LbmControlData *mpControl = this;
00951 #           define TESTGET_FORCE(lev,i,j,k)   mpControl->mCpForces[lev][ ((k*mLevel[lev].lSizey)+j)*mLevel[lev].lSizex+i ]
00952     
00953             glBegin(GL_LINES);
00954             //const int lev=0; 
00955             for(int lev=0; lev<=mMaxRefine; lev++) {
00956             FSGR_FORIJK_BOUNDS(lev) {
00957                     LbmVec pos = LbmVec( 
00958                         ((mvGeoEnd[0]-mvGeoStart[0])/(LbmFloat)mLevel[lev].lSizex) * ((LbmFloat)i+0.5) + mvGeoStart[0], 
00959                         ((mvGeoEnd[1]-mvGeoStart[1])/(LbmFloat)mLevel[lev].lSizey) * ((LbmFloat)j+0.5) + mvGeoStart[1], 
00960                         ((mvGeoEnd[2]-mvGeoStart[2])/(LbmFloat)mLevel[lev].lSizez) * ((LbmFloat)k+0.5) + mvGeoStart[2]  );
00961                     if(LBMDIM==2) pos[2] = ((mvGeoEnd[2]-mvGeoStart[2])*0.5 + mvGeoStart[2]);
00962 
00963                     if((mpControl->mDebugMaxdScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt<=0.) )
00964                     if(TESTGET_FORCE(lev,i,j,k).maxDistance>=0.) 
00965                     if(TESTGET_FORCE(lev,i,j,k).maxDistance<CPF_MAXDINIT ) {
00966                         const float scale = mpControl->mDebugMaxdScale*10001.;
00967                         float dx = TESTGET_FORCE(lev,i,j,k).forceMaxd[0];
00968                         float dy = TESTGET_FORCE(lev,i,j,k).forceMaxd[1];
00969                         float dz = TESTGET_FORCE(lev,i,j,k).forceMaxd[2];
00970                         dx *= scale; dy *= scale; dz *= scale;
00971                         glColor3f( 0,1,0 );
00972                         glVertex3f( pos[0],pos[1],pos[2] );
00973                         glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
00974                     } // */
00975                     if((mpControl->mDebugAttScale>0.) && (TESTGET_FORCE(lev,i,j,k).weightAtt>0.)) {
00976                         const float scale = mpControl->mDebugAttScale*100011.;
00977                         float dx = TESTGET_FORCE(lev,i,j,k).forceAtt[0];
00978                         float dy = TESTGET_FORCE(lev,i,j,k).forceAtt[1];
00979                         float dz = TESTGET_FORCE(lev,i,j,k).forceAtt[2];
00980                         dx *= scale; dy *= scale; dz *= scale;
00981                         glColor3f( 1,0,0 );
00982                         glVertex3f( pos[0],pos[1],pos[2] );
00983                         glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
00984                     } // */
00985                     // why check maxDistance?
00986                     if((mpControl->mDebugVelScale>0.) && (TESTGET_FORCE(lev,i,j,k).maxDistance+TESTGET_FORCE(lev,i,j,k).weightVel>0.)) {
00987                         float scale = mpControl->mDebugVelScale*1.;
00988                         float wvscale = TESTGET_FORCE(lev,i,j,k).weightVel;
00989                         float dx = TESTGET_FORCE(lev,i,j,k).forceVel[0];
00990                         float dy = TESTGET_FORCE(lev,i,j,k).forceVel[1];
00991                         float dz = TESTGET_FORCE(lev,i,j,k).forceVel[2];
00992                         scale *= wvscale;
00993                         dx *= scale; dy *= scale; dz *= scale;
00994                         glColor3f( 0.2,0.2,1 );
00995                         glVertex3f( pos[0],pos[1],pos[2] );
00996                         glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
00997                     } // */
00998                     if((mpControl->mDebugCompavScale>0.) && (TESTGET_FORCE(lev,i,j,k).compAvWeight>0.)) {
00999                         const float scale = mpControl->mDebugCompavScale*1.;
01000                         float dx = TESTGET_FORCE(lev,i,j,k).compAv[0];
01001                         float dy = TESTGET_FORCE(lev,i,j,k).compAv[1];
01002                         float dz = TESTGET_FORCE(lev,i,j,k).compAv[2];
01003                         dx *= scale; dy *= scale; dz *= scale;
01004                         glColor3f( 0.2,0.2,1 );
01005                         glVertex3f( pos[0],pos[1],pos[2] );
01006                         glVertex3f( pos[0]+dx,pos[1]+dy,pos[2]+dz );
01007                     } // */
01008             } // att,maxd
01009             }
01010             glEnd();
01011         }
01012     } // cpssi
01013 
01014     //fprintf(stderr,"BLA\n");
01015     glEnable( GL_LIGHTING ); // dont light lines
01016     glShadeModel(GL_SMOOTH);
01017 }
01018 
01019 #else // LBM_USE_GUI==1
01020 void LbmFsgrSolver::cpDebugDisplay(int dispset) { }
01021 #endif // LBM_USE_GUI==1
01022 
01023