lbm_reference
ParticleSystem.cpp
Go to the documentation of this file.
1 //! \file ParticleSystem.cpp
2 //! \date Mar 12, 2009
3 //! \author Florian Rathgeber
4 
5 #include "../confparser/ConfParser.h"
6 #include "ParticleSystem.h"
7 #include "../lbm/LBM_def.h"
8 
9 namespace particles {
10 
11 ParticleSystem::ParticleSystem ( std::string configFileName )
12  : numParticles_( 0 ),
13  sizeBase_( 0. ),
14  sizeVar_( 0. ),
15  numSprites_( 0 ),
16  device_( 0 ),
17  smgr_( 0 ),
18  drvr_( 0 ) {
19  try {
20 
21  ConfParser p;
22  ConfBlock base = p.parse( configFileName );
23  std::cout << "Parsed configuration file " << configFileName << std::endl;
24 
25  setup( base );
26  solver_.setup( base );
27 
28  srand ( time(NULL) );
29 
30  } catch ( std::exception& e ) {
31  std::cerr << e.what() << std::endl;
32  exit( -1 );
33  } catch ( const char* e ) {
34  std::cerr << e << std::endl;
35  exit( -1 );
36  }
37 }
38 
39 void ParticleSystem::setup( ConfBlock& base ) {
40  try {
41 
42  // Read parameters from config file
43 
44  std::cout << "Setting up ParticleSystem..." << std::endl;
45 
46  // Read domain specification
47  ConfBlock* paramBlock = base.find( "domain" );
48  if ( paramBlock == NULL ) throw "No domain size given.";
49  sizeX_ = paramBlock->getParam<int>( "sizeX" );
50  sizeY_ = paramBlock->getParam<int>( "sizeY" );
51  sizeZ_ = paramBlock->getParam<int>( "sizeZ" );
52  std::cout << "Read domain specification:" << std::endl;
53  std::cout << "sizeX : " << sizeX_ << std::endl;
54  std::cout << "sizeY : " << sizeY_ << std::endl;
55  std::cout << "sizeZ : " << sizeZ_ << std::endl;
56 
57  // Read parameter specification
58  paramBlock = base.find( "parameters" );
59  if ( paramBlock == NULL ) throw "No parameter specification given.";
60  alpha_ = paramBlock->getParam<float>( "alpha" );
61  beta_ = paramBlock->getParam<float>( "beta" );
62  k_ = paramBlock->getParam<float>( "k" );
63  float g_x = paramBlock->getParam<float>( "g_x" );
64  float g_y = paramBlock->getParam<float>( "g_y" );
65  float g_z = paramBlock->getParam<float>( "g_z" );
66  gravity_ = core::vector3df( g_x, g_y, g_z );
67  smokeTemp_ = paramBlock->getParam<float>( "smokeTemp" );
68  ambTemp_ = paramBlock->getParam<float>( "ambTemp" );
69  maxSteps_ = paramBlock->getParam<int>( "maxSteps" );
70  std::cout << "Read parameter specification:" << std::endl;
71  std::cout << "alpha (conservation coefficient) : " << alpha_ << std::endl;
72  std::cout << "beta (transferability coefficient) : " << beta_ << std::endl;
73  std::cout << "k (thermal expansion coefficient) : " << k_ << std::endl;
74  std::cout << "Gravity unit vector : <" << g_x << "," << g_y << "," << g_z << ">" << std::endl;
75  std::cout << "Temperature threshold for smoke (K): " << smokeTemp_ << std::endl;
76  std::cout << "Ambient temperature (K) : " << ambTemp_ << std::endl;
77  std::cout << "Number of steps : " << maxSteps_ << std::endl;
78 
79  // Precompute Gauss function for thermal diffusion
80  int maxElem = sqrt( sizeX_ * sizeX_ + sizeY_ * sizeY_ + sizeZ_ * sizeZ_ );
81  gaussTable_.reserve( maxElem );
82  float a = 1. / sqrt( 2. * M_PI );
83  for ( int i = 0; i < maxElem; ++i ) {
84  gaussTable_.push_back( a * exp(- i * i / 2. ) );
85  }
86 
87  // Set up particle system
88  paramBlock = base.find( "particles" );
89  if ( paramBlock == NULL ) throw "No particle emitters specified.";
90  std::cout << "Set up the particle system..." << std::endl;
91 
92  // Optionally write particle update times to file
93  paramBlock->getParam<std::string>( "updateTimeChart", updFileName_ );
94 
95  ConfBlock::childIterPair cip = paramBlock->findAll( "emitter" );
96 
97  float maxTemp = 0;
98  for ( ConfBlock::childIter it = cip.first; it != cip.second; ++it ) {
99 
100  ConfBlock b = it->second;
101 
102  float xStart = b.getParam<float>( "xStart" );
103  float xEnd = b.getParam<float>( "xEnd" );
104  float yStart = b.getParam<float>( "yStart" );
105  float yEnd = b.getParam<float>( "yEnd" );
106  float zStart = b.getParam<float>( "zStart" );
107  float zEnd = b.getParam<float>( "zEnd" );
108  float temp = b.getParam<float>( "temp" );
109  if ( temp > maxTemp ) maxTemp = temp;
110  int fuel = b.getParam<int>( "fuel" );
111  int emitThreshold = b.getParam<int>( "emitThreshold" );
112  float fuelConsumption = b.getParam<float>( "fuelConsumption" );
113  float lifetimeCoeff = b.getParam<float>( "lifetimeCoeff" );
114  std::cout << "Read emitter specification:" << std::endl;
115  std::cout << "xStart : " << xStart << std::endl;
116  std::cout << "xEnd : " << xEnd << std::endl;
117  std::cout << "yStart : " << yStart << std::endl;
118  std::cout << "yEnd : " << yEnd << std::endl;
119  std::cout << "zStart : " << zStart << std::endl;
120  std::cout << "zEnd : " << zEnd << std::endl;
121  std::cout << "temp : " << temp << std::endl;
122  std::cout << "fuel : " << fuel << std::endl;
123  std::cout << "emitThreshold : " << emitThreshold << std::endl;
124  std::cout << "fuelConsumption : " << fuelConsumption << std::endl;
125  std::cout << "lifetimeCoeff : " << lifetimeCoeff << std::endl;
126 
127  emitters_.push_back( Emitter( core::vector3df( xStart, yStart, zStart ),
128  core::vector3df( xEnd - xStart,
129  yEnd - yStart,
130  zEnd - zStart ),
131  temp,
132  fuel,
133  emitThreshold,
134  fuelConsumption,
135  lifetimeCoeff ) );
136 
137  }
138 
139  // Set up irrlicht engine
140  paramBlock = base.find( "irrlicht" );
141  if ( paramBlock == NULL ) {
142  std::cout << "No irrlicht configuration specified in configuration file. No real-time visualization will be done." << std::endl;
143  } else {
144  std::cout << "Set up irrlicht configuration..." << std::endl;
145 
146  // Read basic visualization parameters
147  sizeBase_ = paramBlock->getParam<float>( "sizeBase" );
148  sizeVar_ = paramBlock->getParam<float>( "sizeVar" );
149  numSprites_ = paramBlock->getParam<int>( "numSprites" );
150  int xRes = paramBlock->getParam<int>( "xRes" );
151  int yRes = paramBlock->getParam<int>( "yRes" );
152  dynamicLights_ = paramBlock->getParam<bool>( "dynamicLights" );
153  std::string camera = paramBlock->getParam<std::string>( "camera" );
154  // Optionally write particle update times to file
155  paramBlock->getParam<std::string>( "irrlichtTimeChart", irrFileName_ );
156  paramBlock->getParam<std::string>( "screenshots", screenshots_ );
157  paramBlock->getParam<int>( "screenStep", screenStep_ );
158 
159  // Create OpenGL device
160  device_ = createDevice( video::EDT_OPENGL, core::dimension2du(xRes,yRes), 24 );
161  if (!device_) throw "Could not create OpenGL device.";
162  smgr_ = device_->getSceneManager();
163  drvr_ = device_->getVideoDriver();
164  assert( smgr_ && drvr_ );
165 
166  // Optionally set ambient light
167  int ambLight;
168  if ( paramBlock->getParam<int>( "ambLight", ambLight ) )
169  smgr_->setAmbientLight( video::SColor(255,ambLight,ambLight,ambLight) );
170 
171  // Add camera to scene
172  if ( camera == "animated" ) {
173  scene::ICameraSceneNode* cam =
174  smgr_->addCameraSceneNode( 0,
175  core::vector3df( sizeX_/2, sizeY_/2, -sizeZ_/2 ),
176  core::vector3df( sizeX_/2, sizeY_/2, 0 ) );
177  cam->addAnimator( smgr_->createFlyCircleAnimator(
178  core::vector3df( sizeX_/2, sizeY_/2, sizeZ_/2 ), // center
179  sizeZ_, // radius
180  0.0002f, // speed
181  core::vector3df(0.f, 1.f, 0.f) // direction
182  ) );
183  } else if ( camera == "fps" ) {
184  scene::ICameraSceneNode* cam =
185  smgr_->addCameraSceneNodeFPS( 0, 100.0f, .1f );
186  cam->setPosition(core::vector3df( sizeX_/2, sizeY_/2, -sizeZ_/2 ));
187  } else {
188  smgr_->addCameraSceneNode( 0,
189  core::vector3df( sizeX_/2, sizeY_/2, -sizeZ_/2 ),
190  core::vector3df( sizeX_/2, sizeY_/2, 0 ) );
191  }
192 
193  cip = paramBlock->findAll( "texture" );
194 
195  for ( ConfBlock::childIter it = cip.first; it != cip.second; ++it ) {
196 
197  ConfBlock b = it->second;
198 
199  std::string file = b.getParam<std::string>( "file" );
200  textures_.push_back( drvr_->getTexture( file.c_str() ) );
201  }
202 
203  cip = paramBlock->findAll( "plane" );
204 
205  for ( ConfBlock::childIter it = cip.first; it != cip.second; ++it ) {
206 
207  ConfBlock b = it->second;
208 
209  std::string name = b.getParam<std::string>( "name" );
210  std::string texture;
211  b.getParam<std::string>( "texture", texture );
212  bool lighting = dynamicLights_;
213  b.getParam<bool>( "lighting", lighting );
214  float sizeTile = b.getParam<float>( "sizeTile" );
215  int numTile = b.getParam<int>( "numTile" );
216  float xCenter = b.getParam<float>( "xCenter" );
217  float yCenter = b.getParam<float>( "yCenter" );
218  float zCenter = b.getParam<float>( "zCenter" );
219  float xRot = 0., yRot = 0., zRot = 0.;
220  b.getParam<float>( "xRot", xRot );
221  b.getParam<float>( "yRot", yRot );
222  b.getParam<float>( "zRot", zRot );
223 
224  // Create a terrain scenenode
225  scene::IAnimatedMesh *terrain_model =
226  smgr_->addHillPlaneMesh( name.c_str(), // Name of the scenenode
227  core::dimension2d<f32>( sizeTile, sizeTile ), // Tile size
228  core::dimension2d<u32>(numTile, numTile), // Tile count
229  0, // Material
230  0.0f, // Hill height
231  core::dimension2d<f32>(0.0f, 0.0f), // countHills
232  core::dimension2d<f32>(numTile, numTile)); // textureRepeatCount
233 
234  scene::IAnimatedMeshSceneNode* terrain_node =
235  smgr_->addAnimatedMeshSceneNode( terrain_model,
236  0,
237  -1,
238  core::vector3df(xCenter,yCenter,zCenter),
239  core::vector3df(xRot,yRot,zRot));
240  if ( texture.length() ) terrain_node->setMaterialTexture(0, drvr_->getTexture( texture.c_str() ));
241  terrain_node->setMaterialFlag(video::EMF_LIGHTING, lighting);
242  }
243 
244  cip = paramBlock->findAll( "mesh" );
245 
246  for ( ConfBlock::childIter it = cip.first; it != cip.second; ++it ) {
247 
248  ConfBlock b = it->second;
249 
250  std::string file = b.getParam<std::string>( "file" );
251  std::string texture0, texture1;
252  b.getParam<std::string>( "texture0", texture0 );
253  b.getParam<std::string>( "texture1", texture1 );
254  bool lighting = dynamicLights_;
255  b.getParam<bool>( "lighting", lighting );
256  float xCenter = b.getParam<float>( "xCenter" );
257  float yCenter = b.getParam<float>( "yCenter" );
258  float zCenter = b.getParam<float>( "zCenter" );
259  float scale = b.getParam<float>( "scale" );
260 
261  scene::IAnimatedMesh* mesh = smgr_->getMesh( file.c_str() );
262 
263  if (mesh) {
264  scene::ISceneNode* node = smgr_->addMeshSceneNode(mesh->getMesh(0));
265  if (node) {
266  node->setPosition(core::vector3df(xCenter,yCenter,zCenter));
267  if ( texture0.length() ) node->setMaterialTexture(0, drvr_->getTexture( texture0.c_str() ));
268  if ( texture1.length() ) node->setMaterialTexture(1, drvr_->getTexture( texture1.c_str() ));
269  node->setMaterialFlag(video::EMF_LIGHTING, lighting);
270  node->setScale( core::vector3df( scale, scale, scale ) );
271  }
272  }
273  }
274 
275  // Generate black body color table
276  std::cout << "Generate black body color table..." << std::endl;
277  generateBlackBodyColorTable( 2. * maxTemp );
278 
279  }
280 
281  // Set up bounding boxes of obstacles
282  paramBlock = base.find( "obstacles" );
283  if ( paramBlock == NULL ) {
284  std::cout << "No obstacles defined." << std::endl;
285  } else {
286 
287  std::cout << "Set up the obstacles..." << std::endl;
288 
289  ConfBlock::childIterPair bit = paramBlock->findAll( "cuboid_stationary" );
290  for ( ConfBlock::childIter it = bit.first; it != bit.second; ++it ) {
291 
292  ConfBlock& b = it->second;
293 
294  bool visible = b.getParam<bool>( "visible" );
295  int xStart = b.getParam<int>( "xStart" );
296  int xEnd = b.getParam<int>( "xEnd" ) + 1;
297  int yStart = b.getParam<int>( "yStart" );
298  int yEnd = b.getParam<int>( "yEnd" ) + 1;
299  int zStart = b.getParam<int>( "zStart" );
300  int zEnd = b.getParam<int>( "zEnd" ) + 1;
301  int dx = xEnd - xStart;
302  int dy = yEnd - yStart;
303  int dz = zEnd - zStart;
304 
305  std::cout << "Stationary cuboid ranging from <" << xStart << ",";
306  std::cout << yStart << "," << zStart << "> to <" << xEnd << "," << yEnd;
307  std::cout << "," << zEnd << ">" << std::endl;
308 
309  if ( visible && device_ ) {
310  std::string texture;
311  b.getParam<std::string>( "texture", texture );
312  bool lighting = dynamicLights_;
313  b.getParam<bool>( "lighting", lighting );
314  scene::IMeshSceneNode* mesh = smgr_->addCubeSceneNode ( 1.0f, // size
315  0, // parent
316  -1, // id
317  core::vector3df( xStart + 0.5 * dx,
318  yStart + 0.5 * dy,
319  zStart + 0.5 * dz ), // position
320  core::vector3df( 0, 0, 0 ), // rotation
321  core::vector3df( dx, dy, dz ) // scale
322  );
323  mesh->setMaterialFlag( video::EMF_LIGHTING, lighting );
324  if ( texture.length() ) mesh->setMaterialTexture( 0, drvr_->getTexture( texture.c_str() ) );
325  }
326  obstacles_.push_back(
327  core::aabbox3df( xStart, yStart, zStart, xEnd, yEnd, zEnd ) );
328  }
329 
330  bit = paramBlock->findAll( "sphere_stationary" );
331  for ( ConfBlock::childIter it = bit.first; it != bit.second; ++it ) {
332 
333  ConfBlock& bl = it->second;
334 
335  bool visible = bl.getParam<bool>( "visible" );
336  float xCenter = bl.getParam<float>( "xCenter" );
337  float yCenter = bl.getParam<float>( "yCenter" );
338  float zCenter = bl.getParam<float>( "zCenter" );
339  float radius = bl.getParam<float>( "radius" );
340 
341  std::cout << "Stationary sphere centered at <" << xCenter << ",";
342  std::cout << yCenter << "," << zCenter << "> with radius " << radius;
343  std::cout << std::endl;
344 
345  // Add scene nodes for moving sphere
346  if ( visible && device_ ) {
347  std::string texture;
348  bl.getParam<std::string>( "texture", texture );
349  bool lighting = dynamicLights_;
350  bl.getParam<bool>( "lighting", lighting );
351  scene::IMeshSceneNode* mesh = smgr_->addSphereSceneNode( radius,
352  128,
353  0,
354  -1,
355  core::vector3df( xCenter, yCenter, zCenter ) );
356  mesh->setMaterialFlag( video::EMF_LIGHTING, lighting );
357  if ( texture.length() ) mesh->setMaterialTexture( 0, drvr_->getTexture( texture.c_str() ) );
358  spheres_.push_back( Sphere( xCenter, yCenter, zCenter, radius, 0.0, 0.0, 0.0, mesh ) );
359  } else {
360  spheres_.push_back( Sphere( xCenter, yCenter, zCenter, radius ) );
361  }
362  }
363 
364  bit = paramBlock->findAll( "sphere_moving" );
365  for ( ConfBlock::childIter it = bit.first; it != bit.second; ++it ) {
366 
367  ConfBlock& bl = it->second;
368 
369  bool visible = bl.getParam<bool>( "visible" );
370  float xCenter = bl.getParam<float>( "xCenter" );
371  float yCenter = bl.getParam<float>( "yCenter" );
372  float zCenter = bl.getParam<float>( "zCenter" );
373  float radius = bl.getParam<float>( "radius" );
374  float u_x = bl.getParam<float>( "u_x" );
375  float u_y = bl.getParam<float>( "u_y" );
376  float u_z = bl.getParam<float>( "u_z" );
377 
378  std::cout << "Moving sphere centered at <" << xCenter << ",";
379  std::cout << yCenter << "," << zCenter << "> with radius " << radius;
380  std::cout << " and u=<" << u_x << "," << u_y << "," << u_z << ">";
381  std::cout << std::endl;
382 
383  // Add scene nodes for moving sphere
384  if ( visible && device_ ) {
385  std::string texture;
386  bl.getParam<std::string>( "texture", texture );
387  bool lighting = dynamicLights_;
388  bl.getParam<bool>( "lighting", lighting );
389  scene::IMeshSceneNode* mesh = smgr_->addSphereSceneNode( radius,
390  128,
391  0,
392  -1,
393  core::vector3df( xCenter, yCenter, zCenter ) );
394  mesh->setMaterialFlag( video::EMF_LIGHTING, lighting );
395  if ( texture.length() ) mesh->setMaterialTexture( 0, drvr_->getTexture( texture.c_str() ) );
396  spheres_.push_back( Sphere( xCenter, yCenter, zCenter, radius, u_x, u_y, u_z, mesh ) );
397  } else {
398  spheres_.push_back( Sphere( xCenter, yCenter, zCenter, radius, u_x, u_y, u_z ) );
399  }
400  }
401  }
402 
403  } catch ( std::exception& e ) {
404  std::cerr << e.what() << std::endl;
405  exit( -1 );
406  } catch ( const char* e ) {
407  std::cerr << e << std::endl;
408  exit( -1 );
409  }
410 
411  std::cout << "ParticleSystem setup finished!" << std::endl;
412 }
413 
414 void ParticleSystem::run() {
415 
416  // Discard first 20 steps to initialize velocity field
417  for ( int i = 0; i < 20; ++i ) solver_.runStep();
418  int step = 20;
419 
420  // Main loop with real-time visualization enabled
421  if ( device_ ) {
422 
423  float totTime = 0., eTimeTot = 0., iTimeTot = 0., sTimeTot = 0., uTimeTot = 0.;
424  std::ofstream updFile, irrFile;
425  if ( updFileName_.length() ) {
426  updFile.open( updFileName_.c_str(), std::ios::out );
427  updFile << "\"Number of Particles\" \"Step time\" \"Emission time\" \"LBM step time\" \"Update time\"\n";
428  }
429  if ( irrFileName_.length() ) {
430  irrFile.open( irrFileName_.c_str(), std::ios::out );
431  irrFile << "\"Number of Particles\" \"Sprites per particle\" \"Step time\" \"Emission time\" \"Rendering time\" \"LBM step time\" \"Update time\" \"FPS\"\n";
432  }
433 
434  // Start the simulation loop
435  while ( device_->run() ) {
436 
437  struct timeval start, end;
438  float stepTime = 0.;
439 
440  gettimeofday(&start, NULL);
441 
442  // emit particles
443  emitParticles();
444 
445  gettimeofday(&end, NULL);
446  float eTime = getTime( start, end );
447 
448  // Draw all primitives
449  drvr_->beginScene(true, true, video::SColor(255,0,0,0));
450  smgr_->drawAll();
451  drvr_->endScene();
452 
453  gettimeofday(&start, NULL);
454  float iTime = getTime( end, start );
455 
456  // Update the window caption
457  core::stringw str = L"LBM fire simulation [";
458  str += drvr_->getName();
459  str += "L], FPS: ";
460  str += (s32)drvr_->getFPS();
461  str += L", Particles: ";
462  str += numParticles_;
463  str += L", Time step: ";
464  str += step;
465  str += L" / ";
466  str += maxSteps_;
467  device_->setWindowCaption( str.c_str() );
468 
469  if ( screenStep_ > 0 && screenshots_.length() && step % screenStep_ == 0 ) {
470  std::ostringstream oss;
471  oss << screenshots_ << "." << std::setw(4) << std::setfill('0') << step << ".jpg";
472  drvr_->writeImageToFile( drvr_->createScreenShot(), oss.str().c_str(), 85 );
473  }
474 
475  // Simulate one LBM step
476  float sTime = solver_.runStep();
477 
478  gettimeofday(&start, NULL);
479 
480  // Update the particles
481  updateParticles();
482 
483  gettimeofday(&end, NULL);
484  float uTime = getTime( start, end );
485 
486  stepTime = eTime + iTime + sTime + uTime;
487  totTime += stepTime;
488  eTimeTot += eTime;
489  iTimeTot += iTime;
490  sTimeTot += sTime;
491  uTimeTot += uTime;
492 
493  std::cout << step << " / " << maxSteps_ << ": ";
494  std::cout << eTime << " + " << iTime << " + " << sTime << " + " << uTime;
495  std::cout << " = " << stepTime << "s, particles: " << numParticles_;
496  std::cout << std::endl;
497 
498  if ( updFileName_.length() )
499  updFile << numParticles_ << " " << stepTime << " " << eTime << " " << sTime << " " << uTime << "\n";
500  if ( irrFileName_.length() )
501  irrFile << numParticles_ << " " << numSprites_ << " " << stepTime << " " << eTime << " " << iTime << " " << sTime << " " << uTime << " " << drvr_->getFPS() << "\n";
502 
503  // Check for end of simulation
504  if ( ++step >= maxSteps_ ) break;
505  }
506 
507  if ( updFileName_.length() ) updFile.close();
508  if ( irrFileName_.length() ) irrFile.close();
509 
510  } else {
511 
512  float totTime = 0.;
513  std::ofstream updFile;
514  if ( updFileName_.length() ) {
515  updFile.open( updFileName_.c_str(), std::ios::out );
516  updFile << "\"Number of Particles\" \"Step time\" \"Emission time\" \"LBM step time\" \"Update time\"\n";
517  }
518  // Start the simulation loop
519  for ( int step = 20; step < maxSteps_; ++step ) {
520 
521  struct timeval start, end;
522  float stepTime = 0.;
523 
524  gettimeofday(&start, NULL);
525 
526  // emit particles
527  emitParticles();
528 
529  gettimeofday(&end, NULL);
530  float eTime = getTime( start, end );
531 
532  // Simulate one LBM step
533  float sTime = solver_.runStep();
534 
535  gettimeofday(&start, NULL);
536 
537  // Update the particles
538  updateParticles();
539 
540  gettimeofday(&end, NULL);
541  float uTime = getTime( start, end );
542 
543  stepTime = eTime + sTime + uTime;
544 
545  std::cout << "Time step " << step << " of " << maxSteps_ << " took ";
546  std::cout << eTime << " + " << sTime << " + " << uTime << " = ";
547  std::cout << stepTime << "secs with " << numParticles_ << " particles";
548  std::cout << std::endl;
549 
550  if ( updFileName_.length() )
551  updFile << numParticles_ << " " << stepTime << " " << eTime << " " << sTime << " " << uTime << "\n";
552  }
553 
554  if ( updFileName_.length() ) updFile.close();
555  }
556 }
557 
558 inline void ParticleSystem::updateParticles() {
559 
560  // Go over all particles of all emitters
561  std::vector< Emitter >::iterator ite, ite2;
562  std::list< Particle >::iterator itp, itp2;
563  for ( ite = emitters_.begin(); ite != emitters_.end(); ++ite ) {
564  // Precompute modificator for particle size
565  float szCoeff = 1. / ( (*ite).fuel_ * (*ite).lifetimeCoeff_ );
566  for ( itp = (*ite).particles_.begin(); itp != (*ite).particles_.end(); ++itp) {
567 
568  // Remove particles that have left the domain or exceeded their lifetime
569  while ( itp != (*ite).particles_.end() && ( (*itp).lifetime_ < 1 ||
570  (*itp).getPos().X < 1 || (*itp).getPos().X > sizeX_ - 1 ||
571  (*itp).getPos().Y < 1 || (*itp).getPos().Y > sizeY_ - 1 ||
572  (*itp).getPos().Z < 1 || (*itp).getPos().Z > sizeZ_ -1 ||
573  (*itp).temp_ < ambTemp_ ) ) {
574 // std::cout << "Particle at pos <" << (*itp).getPos().X << "," << (*itp).getPos().Y << "," << (*itp).getPos().Z <<"> with lifetime " << (*itp).lifetime_ << " removed." << std::endl;
575  (*itp).clear();
576  itp = (*ite).particles_.erase( itp );
577  numParticles_--;
578  }
579 
580  // If no particles are left anymore, go to next emitter
581  if ( itp == (*ite).particles_.end() ) break;
582 
583  // Move particle according to fluid velocity at current position
584  Vec3<float> vel = solver_.getVelocity( (*itp).pos_.X, (*itp).pos_.Y, (*itp).pos_.Z );
585  core::vector3df v = core::vector3df( vel[0], vel[1], vel[2] );
586  // Buoyancy force
587  core::vector3df g = gravity_ * k_ * ( (*itp).temp_ - ambTemp_ );
588  // Check whether the buoyancy force would carry the particle into an
589  // obstacle
590  bool inObstacle = false;
591  core::vector3df u = (*itp).pos_ + g + v;
592  for ( uint i = 0; i < obstacles_.size(); ++i )
593  if( obstacles_[i].isPointInside( u ) ) {
594 // std::cout << "Point <" << u.X << "," << u.Y << "," << u.Z << "> inside obstacle " << i << std::endl;
595  inObstacle = true;
596  break;
597  }
598  for ( uint i = 0; i < spheres_.size(); ++i )
599  if( spheres_[i].isPointInside( u ) ) {
600 // std::cout << "Point <" << u.X << "," << u.Y << "," << u.Z << "> inside obstacle " << i << std::endl;
601  inObstacle = true;
602  break;
603  }
604  // If not, apply it
605  if ( !inObstacle ) (*itp).updatePos( v + g );
606  // If yes, only move it by the fluid velocity
607  else (*itp).updatePos( v );
608 
609  // Update lifetime and particle size
610  float sz = sizeBase_ + sizeVar_ * (*itp).lifetime_ * szCoeff;
611  (*itp).setSize( sz );
612 
613  // Update temperature
614  float tempExt = - gaussTable_[0] * (*itp).temp_;
615  // Add up temperature contributions of all other particles weighted by distance
616  for ( ite2 = emitters_.begin(); ite2 != emitters_.end(); ++ite2 ) {
617  for ( itp2 = (*ite2).particles_.begin(); itp2 != (*ite2).particles_.end(); ++itp2) {
618  tempExt += gaussTable_[ (int) (*itp).dist( *itp2 ) ] * (*itp2).temp_;
619  }
620  }
621  (*itp).temp_ = alpha_ * (*itp).temp_ + beta_ * tempExt;
622 
623  if ( (*itp).type_ == FIRE ) {
624 
625  // If temperature has fallen below threshold, convert to smoke particle
626  if ( (*itp).temp_ < smokeTemp_ ) {
627  (*itp).setSmoke( szCoeff );
628  } else if ( device_ ) {
629  assert( (uint) (((*itp).temp_ - smokeTemp_) / 50.) < bbColorTable_.size() );
630  (*itp).setColor( bbColorTable_[ (int) (((*itp).temp_ - smokeTemp_) / 50.) ] );
631  }
632 
633  } else if ( device_ ){
634  (*itp).setColor( video::SColor( 255 * (*itp).lifetime_ * szCoeff,
635  255 * (*itp).lifetime_ * szCoeff,
636  255 * (*itp).lifetime_ * szCoeff,
637  255 * (*itp).lifetime_ * szCoeff) );
638  }
639 
640  (*itp).lifetime_--;
641  }
642  }
643 
644  // Move the spheres
645  for ( uint i = 0; i < spheres_.size(); ++i ) {
646  spheres_[i].move();
647  }
648 
649 }
650 
651 inline void ParticleSystem::emitParticles() {
652  // Go over all emitters
653  std::vector< Emitter >::iterator ite;
654  for ( ite = emitters_.begin(); ite != emitters_.end(); ++ite ) {
655  // Emit random number of particles up to emit threshold
656  for ( int i = 0; i < (*ite).emitThreshold_; ++i ) {
657  core::vector3df pos = (*ite).pos_ + core::vector3df(
658  ( std::rand() * (*ite).size_.X ) / (float) RAND_MAX,
659  ( std::rand() * (*ite).size_.Y ) / (float) RAND_MAX,
660  ( std::rand() * (*ite).size_.Z ) / (float) RAND_MAX
661  );
662  // If OpenGL output is enabled, also create billboards
663  if ( device_ )
664  (*ite).particles_.push_back(
665  Particle(
666  smgr_,
667  numParticles_,
668  pos,
669  textures_,
670  numSprites_,
671  (*ite).temp_,
672  bbColorTable_[ (int) (((*ite).temp_ - smokeTemp_) / 50.) ],
673  sizeBase_ + sizeVar_,
674  (int) ((*ite).fuel_ * (*ite).lifetimeCoeff_ ),
675  dynamicLights_
676  ) );
677  // Else create only particles without visual representation
678  else
679  (*ite).particles_.push_back(
680  Particle( pos, (*ite).temp_, (int) ((*ite).fuel_ * (*ite).lifetimeCoeff_ ) ) );
681  // Reduce emitter's fuel
682  (*ite).fuel_ *= (*ite).fuelConsumption_;
683  numParticles_++;
684  }
685  }
686 }
687 
688 void ParticleSystem::generateBlackBodyColorTable( float maxTemp ) {
689 
690  float cie_colour_match[81][3] = {
691  {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
692  {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
693  {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
694  {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
695  {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
696  {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
697  {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
698  {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
699  {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
700  {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
701  {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
702  {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
703  {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
704  {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
705  {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
706  {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
707  {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
708  {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
709  {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
710  {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
711  {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
712  {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
713  {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
714  {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
715  {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
716  {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
717  {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
718  };
719 
720  for ( float t = smokeTemp_; t < maxTemp; t += 50. ) {
721 
722  // Calculate x,y,z colors from solar spectrum
723 
724  int i;
725  float lambda, x = 0, y = 0, z = 0, xyz;
726  for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
727  float Me;
728  // Get black body radiation intensity for given temperature and
729  // wavelength
730  float wlm = lambda * 1e-9; // wavelength in meters
731  Me = (3.74183e-16 * pow(wlm, -5.0)) / (exp(1.4388e-2 / (wlm * t)) - 1.0);
732  x += Me * cie_colour_match[i][0];
733  y += Me * cie_colour_match[i][1];
734  z += Me * cie_colour_match[i][2];
735  }
736  xyz = (x + y + z);
737  x /= xyz;
738  y /= xyz;
739  z /= xyz;
740 
741  // Calculate r,g,b colors from x,y,z colors
742 
743  float xr, yr, zr, xg, yg, zg, xb, yb, zb;
744  float xw, yw, zw;
745  float rx, ry, rz, gx, gy, gz, bx, by, bz;
746  float rw, gw, bw;
747 
748  xr = 0.630; yr = 0.340; zr = 1 - (xr + yr);
749  xg = 0.310; yg = 0.595; zg = 1 - (xg + yg);
750  xb = 0.155; yb = 0.070; zb = 1 - (xb + yb);
751 
752  xw = 0.3127; yw = 0.3291; zw = 1 - (xw + yw);
753 
754  // xyz -> rgb matrix, before scaling to white.
755  rx = (yg * zb) - (yb * zg); ry = (xb * zg) - (xg * zb); rz = (xg * yb) - (xb * yg);
756  gx = (yb * zr) - (yr * zb); gy = (xr * zb) - (xb * zr); gz = (xb * yr) - (xr * yb);
757  bx = (yr * zg) - (yg * zr); by = (xg * zr) - (xr * zg); bz = (xr * yg) - (xg * yr);
758 
759  // White scaling factors.
760  // Dividing by yw scales the white luminance to unity, as conventional.
761  rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw;
762  gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw;
763  bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw;
764 
765  // xyz -> rgb matrix, correctly scaled to white.
766  rx = rx / rw; ry = ry / rw; rz = rz / rw;
767  gx = gx / gw; gy = gy / gw; gz = gz / gw;
768  bx = bx / bw; by = by / bw; bz = bz / bw;
769 
770  // rgb of the desired point
771  float r = (rx * x) + (ry * y) + (rz * z);
772  float g = (gx * x) + (gy * y) + (gz * z);
773  float b = (bx * x) + (by * y) + (bz * z);
774 
775  // Constrain to RGB gammut
776  double w;
777  // Amount of white needed is w = - min(0, *r, *g, *b)
778  w = (0 < r) ? 0 : r;
779  w = (w < g) ? w : g;
780  w = (w < b) ? w : b;
781  w = -w;
782  // Add just enough white to make r, g, b all positive.
783  if (w > 0) {
784  r += w; g += w; b += w;
785  }
786 
787  // Normalize
788  float max = (r > b) ? ( (r > g) ? r : g ) : ( (b > g) ? b : g );
789  r *= 255. / max;
790  g *= 255. / max;
791  b *= 255. / max;
792 
793  bbColorTable_.push_back( video::SColor( 25, r, g, b ) );
794  std::cout << "Temperature " << t << "K: RGB <" << r << ", " << g << ", " << b << ">" << std::endl;
795  }
796  std::cout << "Generated black body color table from " << smokeTemp_;
797  std::cout << "K to " << maxTemp << "K (" << bbColorTable_.size() << " values)" << std::endl;
798 }
799 
800 } // namespace particles
Generated on Sun Jan 5 2014 15:28:21 for lbm_reference by   doxygen 1.8.3.1