 |
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