intern/elbeem/intern/simulation_object.cpp

   1 /******************************************************************************
   2  *
   3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
   4  * Copyright 2003-2006 Nils Thuerey
   5  *
   6  * Basic interface for all simulation modules
   7  *
   8  *****************************************************************************/
   9
  10 #include "simulation_object.h"
  11 #include "solver_interface.h"
  12 #include "ntl_bsptree.h"
  13 #include "ntl_ray.h"
  14 #include "ntl_world.h"
  15 #include "solver_interface.h"
  16 #include "particletracer.h"
  17 #include "elbeem.h"
  18
  19 #ifdef _WIN32
  20 #else
  21 #include <sys/time.h>
  22 #endif
  23
  24
  25 //! lbm factory functions
  26 LbmSolverInterface* createSolver();
  27
  28
  29 /******************************************************************************
  30  * Constructor
  31  *****************************************************************************/
  32 SimulationObject::SimulationObject() :
  33         ntlGeometryShader(),
  34         mGeoStart(-100.0), mGeoEnd(100.0),
  35         mpGiTree(NULL), mpGiObjects(NULL),
  36         mpGlob(NULL),
  37         mPanic( false ),
  38         mDebugType( 1 /* =FLUIDDISPNothing*/ ),
  39         mpLbm(NULL), mpParam( NULL ),
  40         mShowSurface(true), mShowParticles(false),
  41         mSelectedCid( NULL ),
  42         mpElbeemSettings( NULL )
  43
  44 {
  45         mpParam = new Parametrizer();
  46         //for(int i=0; i<MAX_DEBDISPSET; i++) { mDebDispSet[i].type  = (i); mDebDispSet[i].on    = false; mDebDispSet[i].scale = 1.0; }
  47
  48         // reset time
  49         mTime                                           = 0.0;
  50 }
  51
  52
  53 /******************************************************************************
  54  * Destructor
  55  *****************************************************************************/
  56 SimulationObject::~SimulationObject()
  57 {
  58         if(mpGiTree)         delete mpGiTree;
  59         if(mpElbeemSettings) delete mpElbeemSettings;
  60         if(mpLbm)            delete mpLbm;
  61         if(mpParam)          delete mpParam;
  62         if(mpParts)          delete mpParts;
  63         debMsgStd("SimulationObject",DM_MSG,"El'Beem Done!\n",10);
  64 }
  65
  66
  67
  68 /*****************************************************************************/
  69 /*! init tree for certain geometry init */
  70 /*****************************************************************************/
  71 void SimulationObject::initGeoTree() {
  72         // unused!! overriden by solver interface
  73         if(mpGlob == NULL) {
  74                 errFatal("SimulationObject::initGeoTree error","Requires globals!", SIMWORLD_INITERROR);
  75                 return;
  76         }
  77         ntlScene *scene = mpGlob->getSimScene();
  78         mpGiObjects = scene->getObjects();
  79
  80         if(mpGiTree != NULL) delete mpGiTree;
  81         char treeFlag = (1<<(mGeoInitId+4));
  82         mpGiTree = new ntlTree( 20, 4, // warning - fixed values for depth & maxtriangles here...
  83                                                                                                 scene, treeFlag );
  84         // unused!! overriden by solver interface
  85 }
  86
  87 /*****************************************************************************/
  88 /*! destroy tree etc. when geometry init done */
  89 /*****************************************************************************/
  90 void SimulationObject::freeGeoTree() {
  91         if(mpGiTree != NULL) delete mpGiTree;
  92 }
  93
  94
  95
  96 // copy & remember settings for later use
  97 void SimulationObject::copyElbeemSettings(elbeemSimulationSettings *settings) {
  98         mpElbeemSettings = new elbeemSimulationSettings;
  99         *mpElbeemSettings = *settings;
 100
 101         mGeoInitId = settings->domainId+1;
 102         debMsgStd("SimulationObject",DM_MSG,"mGeoInitId="<<mGeoInitId<<", domainId="<<settings->domainId, 8);
 103 }
 104
 105 /******************************************************************************
 106  * simluation interface: initialize simulation using the given configuration file
 107  *****************************************************************************/
 108 extern int glob_mpnum;
 109 int SimulationObject::initializeLbmSimulation(ntlRenderGlobals *glob)
 110 {
 111         if(! isSimworldOk() ) return 1;
 112
 113         // already inited?
 114         if(mpLbm) return 0;
 115
 116         mpGlob = glob;
 117         if(!getVisible()) {
 118                 mpAttrs->setAllUsed();
 119                 return 0;
 120         }
 121
 122
 123         mGeoInitId = mpAttrs->readInt("geoinitid", mGeoInitId,"LbmSolverInterface", "mGeoInitId", false);
 124         //mDimension, mSolverType are deprecated
 125         string mSolverType("");
 126         mSolverType = mpAttrs->readString("solver", mSolverType, "SimulationObject","mSolverType", false );
 127
 128         mpLbm = createSolver();
 129   /* check lbm pointer */
 130         if(mpLbm == NULL) {
 131                 errFatal("SimulationObject::initializeLbmSimulation","Unable to init LBM solver! ", SIMWORLD_INITERROR);
 132                 return 2;
 133         }
 134         debMsgStd("SimulationObject::initialized",DM_MSG,"IdStr:"<<mpLbm->getIdString() <<" LBM solver! ", 2);
 135
 136         mpParts = new ParticleTracer();
 137
 138         // for non-param simulations
 139         mpLbm->setParametrizer( mpParam );
 140         mpParam->setAttrList( getAttributeList() );
 141         // not needed.. done in solver_init: mpParam->setSize ... in solver_interface
 142         mpParam->parseAttrList();
 143
 144         mpLbm->setAttrList( getAttributeList() );
 145         mpLbm->setSwsAttrList( getSwsAttributeList() );
 146         mpLbm->parseAttrList();
 147         mpParts->parseAttrList( getAttributeList() );
 148
 149         if(! isSimworldOk() ) return 3;
 150         mpParts->setName( getName() + "_part" );
 151         mpParts->initialize( glob );
 152         if(! isSimworldOk() ) return 4;
 153
 154         // init material settings
 155         string matMc("default");
 156         matMc = mpAttrs->readString("material_surf", matMc, "SimulationObject","matMc", false );
 157         mShowSurface   = mpAttrs->readInt("showsurface", mShowSurface, "SimulationObject","mShowSurface", false );
 158         mShowParticles = mpAttrs->readInt("showparticles", mShowParticles, "SimulationObject","mShowParticles", false );
 159
 160         checkBoundingBox( mGeoStart, mGeoEnd, "SimulationObject::initializeSimulation" );
 161         mpLbm->setLbmInitId( mGeoInitId );
 162         mpLbm->setGeoStart( mGeoStart );
 163         mpLbm->setGeoEnd( mGeoEnd );
 164         mpLbm->setRenderGlobals( mpGlob );
 165         mpLbm->setName( getName() + "_lbm" );
 166         mpLbm->setParticleTracer( mpParts );
 167         if(mpElbeemSettings) {
 168                 // set further settings from API struct init
 169                 if(mpElbeemSettings->outputPath) this->mOutFilename = string(mpElbeemSettings->outputPath);
 170                 mpLbm->initDomainTrafo( mpElbeemSettings->surfaceTrafo );
 171                 mpLbm->setSmoothing(1.0 * mpElbeemSettings->surfaceSmoothing, 1.0 * mpElbeemSettings->surfaceSmoothing);
 172                 mpLbm->setIsoSubdivs(mpElbeemSettings->surfaceSubdivs);
 173                 mpLbm->setSizeX(mpElbeemSettings->resolutionxyz);
 174                 mpLbm->setSizeY(mpElbeemSettings->resolutionxyz);
 175                 mpLbm->setSizeZ(mpElbeemSettings->resolutionxyz);
 176                 mpLbm->setPreviewSize(mpElbeemSettings->previewresxyz);
 177                 mpLbm->setRefinementDesired(mpElbeemSettings->maxRefine);
 178                 mpLbm->setGenerateParticles(mpElbeemSettings->generateParticles);
 179                 // set initial particles
 180                 mpParts->setNumInitialParticles(mpElbeemSettings->numTracerParticles);
 181
 182                 string dinitType = string("no");
 183                 if     (mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_PARTSLIP) dinitType = string("part");
 184                 else if(mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_FREESLIP) dinitType = string("free");
 185                 else /*if(mpElbeemSettings->domainobsType==FLUIDSIM_OBSTACLE_NOSLIP)*/ dinitType = string("no");
 186                 mpLbm->setDomainBound(dinitType);
 187                 mpLbm->setDomainPartSlip(mpElbeemSettings->domainobsPartslip);
 188                 mpLbm->setDumpVelocities(mpElbeemSettings->generateVertexVectors);
 189                 mpLbm->setFarFieldSize(mpElbeemSettings->farFieldSize);
 190                 debMsgStd("SimulationObject::initialize",DM_MSG,"Added domain bound: "<<dinitType<<" ps="<<mpElbeemSettings->domainobsPartslip<<" vv"<<mpElbeemSettings->generateVertexVectors<<","<<mpLbm->getDumpVelocities(), 9 );
 191
 192                 debMsgStd("SimulationObject::initialize",DM_MSG,"Set ElbeemSettings values "<<mpLbm->getGenerateParticles(),10);
 193         }
 194
 195         if(! mpLbm->initializeSolverMemory()   )         { errMsg("SimulationObject::initialize","initializeSolverMemory failed"); mPanic=true; return 10; }
 196         if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverMemory status"); mPanic=true; return 11; }
 197         if(! mpLbm->initializeSolverGrids()    )         { errMsg("SimulationObject::initialize","initializeSolverGrids  failed"); mPanic=true; return 12; }
 198         if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverGrids  status"); mPanic=true; return 13; }
 199         if(! mpLbm->initializeSolverPostinit() )         { errMsg("SimulationObject::initialize","initializeSolverPostin failed"); mPanic=true; return 14; }
 200         if(checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0)) { errMsg("SimulationObject::initialize","initializeSolverPostin status"); mPanic=true; return 15; }
 201
 202         // print cell type stats
 203         bool printStats = true;
 204         if(glob_mpnum>0) printStats=false; // skip in this case
 205         if(printStats) {
 206                 const int jmax = sizeof(CellFlagType)*8;
 207                 int totalCells = 0;
 208                 int flagCount[jmax];
 209                 for(int j=0; j<jmax ; j++) flagCount[j] = 0;
 210                 int diffInits = 0;
 211                 LbmSolverInterface::CellIdentifier cid = mpLbm->getFirstCell();
 212                 for(; mpLbm->noEndCell( cid );
 213                                         mpLbm->advanceCell( cid ) ) {
 214                         int flag = mpLbm->getCellFlag(cid,0);
 215                         int flag2 = mpLbm->getCellFlag(cid,1);
 216                         if(flag != flag2) {
 217                                 diffInits++;
 218                         }
 219                         for(int j=0; j<jmax ; j++) {
 220                                 if( flag&(1<<j) ) flagCount[j]++;
 221                         }
 222                         totalCells++;
 223                 }
 224                 mpLbm->deleteCellIterator( &cid );
 225
 226                 char charNl = '\n';
 227                 debugOutNnl("SimulationObject::initializeLbmSimulation celltype stats: " <<charNl, 5);
 228                 debugOutNnl("no. of cells = "<<totalCells<<", "<<charNl ,5);
 229                 for(int j=0; j<jmax ; j++) {
 230                         std::ostringstream out;
 231                         if(flagCount[j]>0) {
 232                                 out<<"\t" << flagCount[j] <<" x "<< convertCellFlagType2String( (CellFlagType)(1<<j) ) <<", " << charNl;
 233                                 debugOutNnl(out.str(), 5);
 234                         }
 235                 }
 236                 // compute dist. of empty/bnd - fluid - if
 237                 // cfEmpty   = (1<<0), cfBnd  = (1<< 2), cfFluid   = (1<<10), cfInter   = (1<<11),
 238                 if(1){
 239                         std::ostringstream out;
 240                         out.precision(2); out.width(4);
 241                         int totNum = flagCount[1]+flagCount[2]+flagCount[7]+flagCount[8];
 242                         double ebFrac = (double)(flagCount[1]+flagCount[2]) / totNum;
 243                         double flFrac = (double)(flagCount[7]) / totNum;
 244                         double ifFrac = (double)(flagCount[8]) / totNum;
 245                         //???
 246                         out<<"\tFractions: [empty/bnd - fluid - interface - ext. if]  =  [" << ebFrac<<" - " << flFrac<<" - " << ifFrac<<"] "<< charNl;
 247
 248                         if(diffInits > 0) {
 249                                 debMsgStd("SimulationObject::initializeLbmSimulation",DM_MSG,"celltype Warning: Diffinits="<<diffInits<<"!" , 5);
 250                         }
 251                         debugOutNnl(out.str(), 5);
 252                 }
 253         } // cellstats
 254
 255         // might be modified by mpLbm
 256         //mpParts->setStart( mGeoStart );?  mpParts->setEnd( mGeoEnd );?
 257         mpParts->setStart( mpLbm->getGeoStart() );
 258         mpParts->setEnd(   mpLbm->getGeoEnd()   );
 259         mpParts->setCastShadows( false );
 260         mpParts->setReceiveShadows( false );
 261         mpParts->searchMaterial( glob->getMaterials() );
 262
 263         // this has to be inited here - before, the values might be unknown
 264         IsoSurface *surf = mpLbm->getSurfaceGeoObj();
 265         if(surf) {
 266                 surf->setName( "final" ); // final surface mesh
 267                 // warning - this might cause overwriting effects for multiple sims and geom dump...
 268                 surf->setCastShadows( true );
 269                 surf->setReceiveShadows( false );
 270                 surf->searchMaterial( glob->getMaterials() );
 271                 if(mShowSurface) mObjects.push_back( surf );
 272         }
 273
 274 #ifdef ELBEEM_PLUGIN
 275         mShowParticles=1; // for e.g. dumping
 276 #endif // ELBEEM_PLUGIN
 277         if((mpLbm->getGenerateParticles()>0.0)||(mpParts->getNumInitialParticles()>0)) {
 278                 mShowParticles=1;
 279                 mpParts->setDumpParts(true);
 280         }
 281                 //debMsgStd("SimulationObject::init",DM_NOTIFY,"Using envvar ELBEEM_DUMPPARTICLE to set mShowParticles, DEBUG!",1);
 282         //}  // DEBUG ENABLE!!!!!!!!!!
 283         if(mShowParticles) {
 284                 mObjects.push_back(mpParts);
 285         }
 286
 287         // add objects to display for debugging (e.g. levelset particles)
 288         vector<ntlGeometryObject *> debugObjs = mpLbm->getDebugObjects();
 289         for(size_t i=0;i<debugObjs.size(); i++) {
 290                 debugObjs[i]->setCastShadows( false );
 291                 debugObjs[i]->setReceiveShadows( false );
 292                 debugObjs[i]->searchMaterial( glob->getMaterials() );
 293                 mObjects.push_back( debugObjs[i] );
 294                 debMsgStd("SimulationObject::init",DM_NOTIFY,"Added debug obj "<<debugObjs[i]->getName(), 10 );
 295         }
 296         return 0;
 297 }
 298
 299 /*! set current frame */
 300 void SimulationObject::setFrameNum(int num) {
 301         // advance parametrizer
 302         mpParam->setFrameNum(num);
 303 }
 304
 305 /******************************************************************************
 306  * simluation interface: advance simulation another step (whatever delta time that might be)
 307  *****************************************************************************/
 308 void SimulationObject::step( void )
 309 {
 310         if(mpParam->getCurrentAniFrameTime()>0.0) {
 311                 // dont advance for stopped time
 312                 mpLbm->step();
 313                 mTime += mpParam->getTimestep();
 314                 //if(mTime>0.001) { errMsg("DEBUG!!!!!!!!","quit mlsu..."); xit(1); } // PROFILE DEBUG TEST!
 315         }
 316         if(mpLbm->getPanic()) mPanic = true;
 317
 318         checkCallerStatus(FLUIDSIM_CBSTATUS_STEP, 0);
 319         //if((mpElbeemSettings)&&(mpElbeemSettings->runsimCallback)) {
 320                 //int ret = (mpElbeemSettings->runsimCallback)(mpElbeemSettings->runsimUserData, FLUIDSIM_CBSTATUS_STEP, 0);
 321                 //errMsg("runSimulationCallback cbtest1"," "<<this->getName()<<" ret="<<ret);
 322         //}
 323   //debMsgStd("SimulationObject::step",DM_MSG," Sim '"<<mName<<"' stepped to "<<mTime<<" (stept="<<(mpParam->getTimestep())<<", framet="<<getFrameTime()<<") ", 10);
 324 }
 325 /*! prepare visualization of simulation for e.g. raytracing */
 326 void SimulationObject::prepareVisualization( void ) {
 327         if(mPanic) return;
 328         mpLbm->prepareVisualization();
 329 }
 330
 331
 332 /******************************************************************************/
 333 /* get current start simulation time */
 334 double SimulationObject::getStartTime( void ) {
 335         //return mpParam->calculateAniStart();
 336         return mpParam->getAniStart();
 337 }
 338 /* get time for a single animation frame */
 339 double SimulationObject::getFrameTime( int frame ) {
 340         return mpParam->getAniFrameTime(frame);
 341 }
 342 /* get time for a single time step  */
 343 double SimulationObject::getTimestep( void ) {
 344         return mpParam->getTimestep();
 345 }
 346
 347
 348 /******************************************************************************
 349  * return a pointer to the geometry object of this simulation
 350  *****************************************************************************/
 351 //ntlGeometryObject *SimulationObject::getGeometry() { return mpMC; }
 352 vector<ntlGeometryObject *>::iterator
 353 SimulationObject::getObjectsBegin()
 354 {
 355         return mObjects.begin();
 356 }
 357 vector<ntlGeometryObject *>::iterator
 358 SimulationObject::getObjectsEnd()
 359 {
 360         return mObjects.end();
 361 }
 362
 363
 364
 365
 366
 367 /******************************************************************************
 368  * GUI - display debug info
 369  *****************************************************************************/
 370
 371 void SimulationObject::drawDebugDisplay() {
 372 #ifndef NOGUI
 373         if(!getVisible()) return;
 374
 375         //if( mDebugType > (MAX_DEBDISPSET-1) ){ errFatal("SimulationObject::drawDebugDisplay","Invalid debug type!", SIMWORLD_GENERICERROR); return; }
 376         //mDebDispSet[ mDebugType ].on = true;
 377         //errorOut( mDebugType <<"//"<< mDebDispSet[mDebugType].type );
 378         mpLbm->debugDisplay( mDebugType );
 379
 380         //::lbmMarkedCellDisplay<>( mpLbm );
 381         mpLbm->lbmMarkedCellDisplay();
 382 #endif
 383 }
 384
 385 /* GUI - display interactive info  */
 386 void SimulationObject::drawInteractiveDisplay()
 387 {
 388 #ifndef NOGUI
 389         if(!getVisible()) return;
 390         if(mSelectedCid) {
 391                 // in debugDisplayNode if dispset is on is ignored...
 392                 mpLbm->debugDisplayNode( FLUIDDISPGrid, mSelectedCid );
 393         }
 394 #endif
 395 }
 396
 397
 398 /*******************************************************************************/
 399 // GUI - handle mouse movement for selection
 400 /*******************************************************************************/
 401 void SimulationObject::setMousePos(int x,int y, ntlVec3Gfx org, ntlVec3Gfx dir)
 402 {
 403         normalize( dir );
 404         // assume 2D sim is in XY plane...
 405
 406         double zplane = (mGeoEnd[2]-mGeoStart[2])*0.5;
 407         double zt = (zplane-org[2]) / dir[2];
 408         ntlVec3Gfx pos(
 409                         org[0]+ dir[0] * zt,
 410                         org[1]+ dir[1] * zt, 0.0);
 411
 412         mSelectedCid = mpLbm->getCellAt( pos );
 413         //errMsg("SMP ", mName<< x<<" "<<y<<" - "<<dir );
 414         x = y = 0; // remove warning
 415 }
 416
 417
 418 void SimulationObject::setMouseClick()
 419 {
 420         if(mSelectedCid) {
 421                 //::debugPrintNodeInfo<>( mpLbm, mSelectedCid, mpLbm->getNodeInfoString() );
 422                 mpLbm->debugPrintNodeInfo( mSelectedCid );
 423         }
 424 }
 425
 426 /*! notify object that dump is in progress (e.g. for field dump) */
 427 void SimulationObject::notifyShaderOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
 428         if(!mpLbm) return;
 429
 430         mpLbm->notifySolverOfDump(dumptype, frameNr,frameNrStr,outfilename);
 431         checkCallerStatus(FLUIDSIM_CBSTATUS_NEWFRAME, frameNr);
 432 }
 433
 434 /*! check status (e.g. stop/abort) from calling program, returns !=0 if sth. happened... */
 435 int SimulationObject::checkCallerStatus(int status, int frame) {
 436         //return 0; // DEBUG
 437         int ret = 0;
 438         if((mpElbeemSettings)&&(mpElbeemSettings->runsimCallback)) {
 439                 ret = (mpElbeemSettings->runsimCallback)(mpElbeemSettings->runsimUserData, status,frame);
 440                 if(ret!=FLUIDSIM_CBRET_CONTINUE) {
 441                         if(ret==FLUIDSIM_CBRET_STOP) {
 442                                 debMsgStd("SimulationObject::notifySolverOfDump",DM_NOTIFY,"Got stop signal from caller",1);
 443                                 setElbeemState( SIMWORLD_STOP );
 444                         }
 445                         else if(ret==FLUIDSIM_CBRET_ABORT) {
 446                                 errFatal("SimulationObject::notifySolverOfDump","Got abort signal from caller, aborting...", SIMWORLD_GENERICERROR );
 447                                 mPanic = 1;
 448                         }
 449                         else {
 450                                 errMsg("SimulationObject::notifySolverOfDump","Invalid callback return value: "<<ret<<", ignoring... ");
 451                         }
 452                 }
 453         }
 454
 455         //debMsgStd("SimulationObject::checkCallerStatus",DM_MSG, "s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
 456         if(isSimworldOk()) return 0;
 457         return 1;
 458 }
 459