.. _program_listing_file_SPlisHSPlasH_Emitter.cpp:

Program Listing for File Emitter.cpp
====================================

|exhale_lsh| :ref:`Return to documentation for file <file_SPlisHSPlasH_Emitter.cpp>` (``SPlisHSPlasH/Emitter.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "Emitter.h"
   #include "SPHKernels.h"
   #include <iostream>
   #include "TimeManager.h"
   #include "TimeStep.h"
   #include "FluidModel.h"
   #include "Simulation.h"
   
   using namespace SPH;
   
   Emitter::Emitter(FluidModel *model,
       const unsigned int width, const unsigned int height,
       const Vector3r &pos, const Matrix3r & rotation,
       const Real velocity,
       const unsigned int type)
       : m_model(model)
       , m_width(width)
       , m_height(height)
       , m_x(pos)
       , m_rotation(rotation)
       , m_velocity(velocity)
       , m_type(type)
       , m_nextEmitTime(0)
       , m_emitStartTime(0)
       , m_emitEndTime(std::numeric_limits<Real>::max())
       , m_emitCounter(0)
   {
   }
   
   Emitter::~Emitter(void)
   {
   }
   
   void Emitter::reset()
   {
       m_nextEmitTime = m_emitStartTime;
       m_emitCounter = 0;  
   }
   
   Vector3r Emitter::getSize(const Real width, const Real height, const int type)
   {
       Simulation *sim = Simulation::getCurrent();
       const Real radius = sim->getParticleRadius();
       const Real diam = static_cast<Real>(2.0)*radius;
   
       const Real supportRadius = sim->getSupportRadius();
       const Real animationMarginAround = diam;
       Vector3r size;
       if (type == 0)
       {
           if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
           {
               size = {
                   2 * supportRadius,
                   height * diam + 2 * animationMarginAround,
                   width * diam + 2 * animationMarginAround
               };
           }
           else
           {
               size = {
                   static_cast<Real>(2.0)* supportRadius,
                   height * diam + static_cast<Real>(2.5) * animationMarginAround,
                   width * diam + static_cast<Real>(2.5) * animationMarginAround
               };
           }
       }
       else
       {
           if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
           {
               // height and radius of cylinder
               const Real h = 2 * supportRadius;
               const Real r = 0.5f * width * diam + animationMarginAround;
               size = { h, 2 * r, 2 * r };
           }
           else
           {
               // height and radius of cylinder
               const Real h = static_cast<Real>(2.25) * supportRadius;
               const Real r = 0.5f * width * diam + animationMarginAround;
               size = { h, static_cast<Real>(2.25) * r, static_cast<Real>(2.25) * r };
           }
       }
   
       return size;
   }
   
   void Emitter::emitParticles(std::vector <unsigned int> &reusedParticles, unsigned int &indexReuse, unsigned int &numEmittedParticles)
   {
       TimeManager *tm = TimeManager::getCurrent();
       const Real t = tm->getTime();
       const Real timeStepSize = tm->getTimeStepSize();
       const Vector3r & emitDir = m_rotation.col(0);
       Vector3r emitVel = m_velocity * emitDir;
       Simulation *sim = Simulation::getCurrent();
       const Real radius = sim->getParticleRadius();
       const Real diam = static_cast<Real>(2.0)*radius;
   
       // shortly before the emitter starts, cleanup the emitter from particles
       if (t < m_emitStartTime || t > m_emitEndTime)
           emitVel = emitDir * radius * 10 / 0.25;
   
       if (t >= m_emitStartTime - 0.25 && t <= m_emitEndTime)
       {
           // animate emitted particles
           const Vector3r & x0 = m_x;
   
           const Real animationMarginAhead = sim->getSupportRadius();
           const Vector3r size = getSize(static_cast<Real>(m_width), static_cast<Real>(m_height), m_type);
           const Vector3r halfSize = 0.5 * size;
           const Vector3r pos = x0 + static_cast<Real>(0.5) * animationMarginAhead * emitDir;
   
           const unsigned int nModels = sim->numberOfFluidModels();
           for (unsigned int m = 0; m < nModels; m++)
           {
               FluidModel *fm = sim->getFluidModel(m);
               const unsigned int numParticles = fm->numActiveParticles();
               #pragma omp parallel for schedule(static) default(shared)
               for (int i = 0; i < (int)numParticles; i++)
               {
                   Vector3r &xi = fm->getPosition(i);
                   if (inBox(xi, pos, m_rotation, halfSize))
                   {
                       fm->getVelocity(i) = emitVel;
                       fm->getPosition(i) += timeStepSize * emitVel;
                       fm->setParticleState(i, ParticleState::AnimatedByEmitter);
                       fm->setObjectId(i, m_objectId);
                   }
               }
           }
       }
       if (t < m_nextEmitTime || t > m_emitEndTime)
       {
           return;
       }
   
       const Vector3r axisHeight = m_rotation.col(1);
       const Vector3r axisWidth = m_rotation.col(2);
       
       const Real startX = -static_cast<Real>(0.5)*(m_width  - 1)*diam;
       const Real startZ = -static_cast<Real>(0.5)*(m_height - 1)*diam;
   
       // t-m_nextEmitTime is the time that has passed between the time the particle has been emitted and the current time step.
       // timeStepSize is added because emission happens at the end of the time step, but the particles are not animated anymore.
       const Real dt = t - m_nextEmitTime + timeStepSize;
   
       const Vector3r velocityOffset = dt * emitVel;
       const Vector3r offset = m_x + velocityOffset;
   
       if ((m_model->numActiveParticles() < m_model->numParticles()) ||
           (reusedParticles.size() > 0))
       {
           unsigned int indexNotReuse = m_model->numActiveParticles();
           for (unsigned int i = 0; i < m_width; i++)
           {
               for (unsigned int j = 0; j < m_height; j++)
               {
                   unsigned int index = 0;
                   bool reused = false;
                   if (indexReuse < reusedParticles.size())
                   {
                       index = reusedParticles[indexReuse];
                       reused = true;
                   }
                   else
                   {
                       index = indexNotReuse;
                   }
   
                   if (index < m_model->numParticles())
                   {
                       m_model->getPosition(index) = (i*diam + startX)*axisWidth + (j*diam + startZ)*axisHeight + offset;
                       m_model->getVelocity(index) = emitVel;
                       m_model->setParticleState(index, ParticleState::AnimatedByEmitter);
                       m_model->setObjectId(index, m_objectId);
   
                       if (reused)
                       {
                           indexReuse++;
                       }
                       else
                       {
                           numEmittedParticles++;
                           indexNotReuse++;
                       }
                       index++;
                   }
               }
           }
   
           if (numEmittedParticles != 0)
           {
               m_model->setNumActiveParticles(m_model->numActiveParticles() + numEmittedParticles);
               sim->emittedParticles(m_model, m_model->numActiveParticles() - numEmittedParticles);
               sim->getNeighborhoodSearch()->resize_point_set(m_model->getPointSetIndex(), &m_model->getPosition(0)[0], m_model->numActiveParticles());
           }
       }
       else
       {
           if (m_model->numActiveParticles() < m_model->numParticles())
           {
               unsigned int index = m_model->numActiveParticles();
               for (unsigned int i = 0; i < m_width; i++)
               {
                   for (unsigned int j = 0; j < m_height; j++)
                   {
                       if (index < m_model->numParticles())
                       {
                           m_model->getPosition(index) = (i*diam + startX)*axisWidth + (j*diam + startZ)*axisHeight + offset;
                           m_model->getVelocity(index) = emitVel;
                           m_model->setParticleState(index, ParticleState::AnimatedByEmitter);
                           m_model->setObjectId(index, m_objectId);
                           numEmittedParticles++;
                       }
                       index++;
                   }
               }
               m_model->setNumActiveParticles(m_model->numActiveParticles() + numEmittedParticles);
               sim->emittedParticles(m_model, m_model->numActiveParticles() - numEmittedParticles);
               sim->getNeighborhoodSearch()->resize_point_set(m_model->getPointSetIndex(), &m_model->getPosition(0)[0], m_model->numActiveParticles());
           }
       }
   
       m_nextEmitTime += diam / m_velocity;
       m_emitCounter++;
   }
   
   
   
   void Emitter::emitParticlesCircle(std::vector <unsigned int> &reusedParticles, unsigned int &indexReuse, unsigned int &numEmittedParticles)
   {
       TimeManager *tm = TimeManager::getCurrent();
       const Real t = tm->getTime();
       const Real timeStepSize = tm->getTimeStepSize();
       const Vector3r & emitDir = m_rotation.col(0);
       Simulation *sim = Simulation::getCurrent();
       const Real particleRadius = sim->getParticleRadius();
       const Real diam = static_cast<Real>(2.0)*particleRadius;
       Vector3r emitVel = m_velocity * emitDir;
   
       // shortly before the emitter starts, cleanup the emitter from particles
       if (t < m_emitStartTime || t > m_emitEndTime)
           emitVel = emitDir * particleRadius * 10 / 0.25;
   
       if (t >= m_emitStartTime - 0.25 && t <= m_emitEndTime)
       {
           // animate emitted particles        
           const Vector3r & x0 = m_x;
   
           const Real animationMarginAhead = sim->getSupportRadius();
           const Vector3r size = getSize(static_cast<Real>(m_width), static_cast<Real>(m_height), m_type);
           const Real h = size[0];
           const Real r = 0.5f * size[1];
           const Real r2 = r * r;
           const Vector3r pos = x0 + 0.5f * animationMarginAhead * emitDir;
   
   
           const unsigned int nModels = sim->numberOfFluidModels();
           for (unsigned int m = 0; m < nModels; m++)
           {
               FluidModel *fm = sim->getFluidModel(m);
               const unsigned int numParticles = fm->numActiveParticles();
               #pragma omp parallel for schedule(static) default(shared)
               for (int i = 0; i < (int)numParticles; i++)
               {
                   Vector3r &xi = fm->getPosition(i);
                   if (inCylinder(xi, pos, m_rotation, h, r2))
                   {
                       fm->getVelocity(i) = emitVel;
                       fm->getPosition(i) += timeStepSize * emitVel;
                       fm->setParticleState(i, ParticleState::AnimatedByEmitter);
                       fm->setObjectId(i, m_objectId);
                   }
               }
           }
       }
       if (t < m_nextEmitTime || t > m_emitEndTime)
       {
           return;
       }
   
       Vector3r axisHeight = m_rotation.col(1);
       Vector3r axisWidth = m_rotation.col(2);
   
       const Real radius = (static_cast<Real>(0.5) * (Real)m_width * diam);
       const Real radius2 = radius*radius;
   
       const Real startX = -static_cast<Real>(0.5)*(m_width - 1)*diam;
       const Real startZ = -static_cast<Real>(0.5)*(m_width - 1)*diam;
   
       // t-m_nextEmitTime is the time that has passed between the time the particle has been emitted and the current time step.
       // timeStepSize is added because emission happens at the end of the time step, but the particles are not animated anymore.
       const Real dt = t - m_nextEmitTime + timeStepSize;
       const Vector3r velocity = emitVel;
       const Vector3r velocityOffset = dt * velocity;
       const Vector3r offset = m_x + velocityOffset;
   
       if ((m_model->numActiveParticles() < m_model->numParticles()) ||
           (reusedParticles.size() > 0))
       {
           unsigned int indexNotReuse = m_model->numActiveParticles();
           for (unsigned int i = 0; i < m_width; i++)
           {
               for (unsigned int j = 0; j < m_width; j++)
               {
                   const Real x = (i*diam + startX);
                   const Real y = (j*diam + startZ);
   
                   unsigned int index = 0;
                   bool reused = false;
                   if (indexReuse < reusedParticles.size())
                   {
                       index = reusedParticles[indexReuse];
                       reused = true;
                   }
                   else
                   {
                       index = indexNotReuse;
                   }
   
                   if ((index < m_model->numParticles()) && (x*x + y*y <= radius2))
                   {
                       m_model->getPosition(index) = x*axisWidth + y*axisHeight + offset;
                       m_model->getVelocity(index) = velocity;
                       m_model->setParticleState(index, ParticleState::AnimatedByEmitter);
                       m_model->setObjectId(index, m_objectId);
   
                       if (reused)
                       {
                           indexReuse++;
                       }
                       else
                       {
                           numEmittedParticles++;
                           indexNotReuse++;
                       }
                       index++;
                   }
               }
           }
   
           if (numEmittedParticles != 0)
           {
               m_model->setNumActiveParticles(m_model->numActiveParticles() + numEmittedParticles);
               sim->emittedParticles(m_model, m_model->numActiveParticles() - numEmittedParticles);
               sim->getNeighborhoodSearch()->resize_point_set(m_model->getPointSetIndex(), &m_model->getPosition(0)[0], m_model->numActiveParticles());
           }
       }
       else
       {
           if (m_model->numActiveParticles() < m_model->numParticles())
           {
               unsigned int index = m_model->numActiveParticles();
               for (unsigned int i = 0; i < m_width; i++)
               {
                   for (unsigned int j = 0; j < m_width; j++)
                   {
                       const Real x = (i*diam + startX);
                       const Real y = (j*diam + startZ);
                       if ((index < m_model->numParticles()) && (x*x + y*y <= radius2))
                       {
                           m_model->getPosition(index) = x*axisWidth + y*axisHeight + offset;
                           m_model->getVelocity(index) = velocity;
                           m_model->setParticleState(index, ParticleState::AnimatedByEmitter);
                           m_model->setObjectId(index, m_objectId);
                           numEmittedParticles++;
                           index++;
                       }
                   }
               }
               m_model->setNumActiveParticles(m_model->numActiveParticles() + numEmittedParticles);
               sim->emittedParticles(m_model, m_model->numActiveParticles() - numEmittedParticles);
               sim->getNeighborhoodSearch()->resize_point_set(m_model->getPointSetIndex(), &m_model->getPosition(0)[0], m_model->numActiveParticles());
           }
       }
   
       m_nextEmitTime += diam / m_velocity;
       m_emitCounter++;
   }
   
   void Emitter::step(std::vector <unsigned int> &reusedParticles, unsigned int &indexReuse, unsigned int &numEmittedParticles)
   {
       if (m_type == 1)
           emitParticlesCircle(reusedParticles, indexReuse, numEmittedParticles);
       else
           emitParticles(reusedParticles, indexReuse, numEmittedParticles);
   }
   
   void SPH::Emitter::saveState(BinaryFileWriter &binWriter)
   {
       binWriter.write(m_nextEmitTime);
       binWriter.write(m_emitCounter);
   }
   
   void SPH::Emitter::loadState(BinaryFileReader &binReader)
   {
       binReader.read(m_nextEmitTime);
       binReader.read(m_emitCounter);
   }