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