#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);

}