https://github.com/InteractiveComputerGraphics/SPlisHSPlasH
Tip revision: 619a911372dca5f48c562d39c13b9f83eda6e10a authored by Jan Bender on 19 January 2022, 09:38:43 UTC
- fixed readme
- fixed readme
Tip revision: 619a911
Emitter.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);
}