Program Listing for File TimeStepPF.cpp

↰ Return to documentation for file (SPlisHSPlasH/PF/TimeStepPF.cpp)

#include "TimeStepPF.h"

#include "SimulationDataPF.h"
#include "SPlisHSPlasH/SPHKernels.h"
#include "SPlisHSPlasH/TimeManager.h"
#include "Utilities/Timing.h"
#include "SPlisHSPlasH/Simulation.h"
#include "SPlisHSPlasH/BoundaryModel_Akinci2012.h"
#include "SPlisHSPlasH/BoundaryModel_Koschier2017.h"
#include "SPlisHSPlasH/BoundaryModel_Bender2019.h"

#include <atomic>
#include <iostream>

using namespace SPH;
using namespace std;
using namespace GenParam;

std::string TimeStepPF::METHOD_NAME = "Projective Fluids";
int TimeStepPF::SOLVER_ITERATIONS = -1;
int TimeStepPF::MIN_ITERATIONS = -1;
int TimeStepPF::MAX_ITERATIONS = -1;
int TimeStepPF::MAX_ERROR = -1;
int TimeStepPF::STIFFNESS = -1;

#define Vec3Block(i) template segment<3>(3 * (i))

// helper functions
namespace
{
    template <typename T>
    struct atomic_wrapper
    {
        std::atomic<T> _a;

        atomic_wrapper() : _a(0) {}
        atomic_wrapper(const std::atomic<T> &a) :_a(a.load()) {}
        atomic_wrapper(const atomic_wrapper &other) :_a(other._a.load()) {}
        atomic_wrapper &operator=(const atomic_wrapper &other)
        {
            _a.store(other._a.load());
            return *this;
        }
    };

    inline void addToAtomicReal(std::atomic<Real> & a, const Real & r)
    {
        Real current = a;
        while (!a.compare_exchange_weak(current, current + r))
            ;
    }
}

using AtomicRealVec = std::vector < atomic_wrapper<Real> >;

TimeStepPF::TimeStepPF() :
    TimeStep(),
    m_stiffness(50000.0),
    m_numActiveParticlesTotal(0)
{
    m_simulationData.init();

    m_iterations = 0;
    m_minIterations = 2;
    m_maxIterations = 100;
    m_maxError = static_cast<Real>(1e-10);

    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        model->addField({ "oldPosition", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getOldPosition(fluidModelIndex, i)[0]; }, true });
        model->addField({ "S", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getS(fluidModelIndex, i)[0]; } });
        model->addField({ "numFluidNeighbors", METHOD_NAME, FieldType::UInt, [this, fluidModelIndex](const unsigned int i) -> unsigned int* { return &m_simulationData.getNumFluidNeighbors(fluidModelIndex, i); } });
        model->addField({ "diag", METHOD_NAME, FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDiag(fluidModelIndex, i)[0]; } });
    }
}

TimeStepPF::~TimeStepPF(void)
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        model->removeFieldByName("oldPosition");
        model->removeFieldByName("S");
        model->removeFieldByName("numFluidNeighbors");
        model->removeFieldByName("diag");
    }
}

void TimeStepPF::initParameters()
{
    TimeStep::initParameters();

    SOLVER_ITERATIONS = createNumericParameter("iterations", "Iterations", &m_iterations);
    setGroup(SOLVER_ITERATIONS, "Simulation|Projective Fluids");
    setDescription(SOLVER_ITERATIONS, "Iterations required by the pressure solver.");
    getParameter(SOLVER_ITERATIONS)->setReadOnly(true);

    MIN_ITERATIONS = createNumericParameter("minIterations", "Min. iterations", &m_minIterations);
    setGroup(MIN_ITERATIONS, "Simulation|Projective Fluids");
    setDescription(MIN_ITERATIONS, "Minimal number of iterations of the pressure solver.");
    static_cast<NumericParameter<unsigned int>*>(getParameter(MIN_ITERATIONS))->setMinValue(0);

    MAX_ITERATIONS = createNumericParameter("maxIterations", "Max. iterations", &m_maxIterations);
    setGroup(MAX_ITERATIONS, "Simulation|Projective Fluids");
    setDescription(MAX_ITERATIONS, "Maximal number of iterations of the pressure solver.");
    static_cast<NumericParameter<unsigned int>*>(getParameter(MAX_ITERATIONS))->setMinValue(1);

    MAX_ERROR = createNumericParameter("maxError", "Max. error", &m_maxError);
    setGroup(MAX_ERROR, "Simulation|Projective Fluids");
    setDescription(MAX_ERROR, "Maximal error.");
    static_cast<RealParameter*>(getParameter(MAX_ERROR))->setMinValue(static_cast<Real>(1e-12));

    STIFFNESS = createNumericParameter("stiffness", "Stiffness", &m_stiffness);
    setGroup(STIFFNESS, "Simulation|Projective Fluids");
    setDescription(STIFFNESS, "Stiffness coefficient.");
    static_cast<RealParameter*>(getParameter(STIFFNESS))->setMinValue(static_cast<Real>(1e-6));
}

void TimeStepPF::step()
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    TimeManager *tm = TimeManager::getCurrent ();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        clearAccelerations(fluidModelIndex);
        initialGuessForPositions(fluidModelIndex);
    }
    sim->performNeighborhoodSearch();

#ifdef USE_PERFORMANCE_OPTIMIZATION
    precomputeValues();
#endif

    if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
        computeVolumeAndBoundaryX();
    else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
        computeDensityAndGradient();

    START_TIMING("solvePDConstraints");
    solvePDConstraints();
    STOP_TIMING_AVG;

    if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
        computeVolumeAndBoundaryX();
    else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
        computeDensityAndGradient();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
        computeDensities(fluidModelIndex);
    sim->computeNonPressureForces();
    addAccellerationToVelocity();

    // update emitters
    sim->emitParticles();
    sim->animateParticles();
    // Compute new time
    sim->updateTimeStepSize();
    tm->setTime(tm->getTime () + tm->getTimeStepSize());
}

void TimeStepPF::reset()
{
    TimeStep::reset();
    m_simulationData.reset();
    m_iterations = 0;
}

void TimeStepPF::initialGuessForPositions(const unsigned int fluidModelIndex)
{
    Simulation *sim = Simulation::getCurrent();
    FluidModel *model = sim->getFluidModel(fluidModelIndex);
    const auto numParticles = model->numActiveParticles();
    const auto h = TimeManager::getCurrent()->getTimeStepSize();

    #pragma omp parallel for
    for (int i = 0; i < (int)numParticles; i++)
    {
        m_simulationData.setOldPosition(fluidModelIndex, i, model->getPosition(i));
        if (model->getParticleState(i) == ParticleState::Active)
        {
            const auto newPos = (model->getPosition(i) + h * model->getVelocity(i) + (h * h) * model->getAcceleration(i)).eval();
            model->setPosition(i, newPos);
            m_simulationData.setS(fluidModelIndex, i, newPos);
        }
        else
        {
            m_simulationData.setS(fluidModelIndex, i, model->getPosition(i));
        }
    }
}

void TimeStepPF::solvePDConstraints()
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nFluids = sim->numberOfFluidModels();
    if (nFluids == 0)
        return;

    // total number of active fluid particles
    m_numActiveParticlesTotal = m_simulationData.getParticleOffset(nFluids - 1) + sim->getFluidModel(nFluids - 1)->numActiveParticles();

    VectorXr x(3 * m_numActiveParticlesTotal);
    VectorXr b(3 * m_numActiveParticlesTotal);

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const unsigned int offset = m_simulationData.getParticleOffset(fluidModelIndex);

        #pragma omp parallel for schedule(static)
        for (int i = 0; i < (int)model->numActiveParticles(); i++)
        {
            // initialize positions
            x.Vec3Block(offset + i) = m_simulationData.getS(fluidModelIndex, i);

            // count number of fluid neighbors for relaxation
            unsigned int nNeighbors = 0;
            for (unsigned int pid = 0; pid < nFluids; pid++)
                nNeighbors += sim->numberOfNeighbors(fluidModelIndex, pid, i);
            m_simulationData.setNumFluidNeighbors(fluidModelIndex, i, nNeighbors + 1u);
        }
    }

    // Init linear system solver and preconditioner
    MatrixReplacement A(3 * m_numActiveParticlesTotal, matrixVecProd, (void*) this);
#ifdef PD_USE_DIAGONAL_PRECONDITIONER
    preparePreconditioner();
    m_solver.preconditioner().init(m_numActiveParticlesTotal, diagonalMatrixElement, (void*)this);
#endif
    m_solver.setMaxIterations(m_maxIterations);
    m_solver.compute(A);

    for (m_iterations = 0u; m_iterations < m_maxIterations; m_iterations++)
    {
        // Compute RHS
        matrixFreeRHS(x,b);

        // Solve linear system
#ifdef PD_USE_DIAGONAL_PRECONDITIONER
        // hack to make the solver perform at least min_iter CG iterations for stability
        const unsigned int min_iter = m_minIterations;
        m_solver.setTolerance(m_iterations < min_iter ? static_cast<Real>(1e-32) : static_cast<Real>(1e-10));
#endif
        x = m_solver.solveWithGuess(b, x);
        if (m_solver.iterations() == 0)
            break;
    }

    updatePositionsAndVelocity(x);
}

void SPH::TimeStepPF::updatePositionsAndVelocity(const VectorXr & x)
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    const Real h_inv = 1 / TimeManager::getCurrent()->getTimeStepSize();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const unsigned int numParticles = model->numActiveParticles();
        const unsigned int offset = m_simulationData.getParticleOffset(fluidModelIndex);

        #pragma omp parallel default(shared)
        {
            #pragma omp for schedule(static)
            for (int i = 0; i < (int)numParticles; i++)
            {
                if (model->getParticleState(i) == ParticleState::Active)
                {
                    const Vector3r vel = h_inv * (x.Vec3Block(offset + i) - m_simulationData.getOldPosition(fluidModelIndex, i));
                    model->setPosition(i, x.Vec3Block(offset + i));
                    model->setVelocity(i, vel);
                }
            }
        }
    }
}

#ifdef PD_USE_DIAGONAL_PRECONDITIONER
FORCE_INLINE void SPH::TimeStepPF::diagonalMatrixElement(const unsigned int row, Vector3r &result, void *userData)
{
    TimeStepPF *timeStepPF = static_cast<TimeStepPF*>(userData);
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();

    // TODO: use m_simulationData->getParticleOffset instead of accessing FluidModels
    // find corresponding fluid model and particle index for the current row
    unsigned int fluidRow = row;
    unsigned int fluidModelIndex;
    for (fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        if (fluidRow >= model->numActiveParticles())
            fluidRow -= model->numActiveParticles();
        else
            break;
    }

    result = timeStepPF->m_simulationData.getDiag(fluidModelIndex, fluidRow);
}

void SPH::TimeStepPF::preparePreconditioner()
{
    Simulation *sim = Simulation::getCurrent();
    const unsigned int nModels = sim->numberOfFluidModels();
    const auto h = TimeManager::getCurrent()->getTimeStepSize();
    const auto system_scale = h * h * m_stiffness;

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const auto numParticles = model->numActiveParticles();
        const Real density0_inv = 1 / model->getDensity0();

        #pragma omp parallel for schedule(static)
        for (int i = 0; i < (int)numParticles; i++)
        {
            Real density_scale = 1;
            for (unsigned int pid = 0; pid < nModels; pid++)
            {
                density_scale += density0_inv * sim->getFluidModelFromPointSet(pid)->getDensity0() * sim->numberOfNeighbors(fluidModelIndex, pid, i);
            }
            m_simulationData.setDiag(fluidModelIndex, i, Vector3r::Constant(system_scale * density_scale + model->getMass(i)));
        }
    }
}

#endif

void SPH::TimeStepPF::addAccellerationToVelocity()
{
    Simulation *sim = Simulation::getCurrent();
    const auto  h = TimeManager::getCurrent()->getTimeStepSize();
    const unsigned int nModels = sim->numberOfFluidModels();

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const auto  numParticles = model->numActiveParticles();
        #pragma omp parallel default(shared)
        {
            #pragma omp for schedule(static)
            for (int i = 0; i < (int)numParticles; i++)
            {
                if (model->getParticleState(i) == ParticleState::Active)
                    model->setVelocity(i, model->getVelocity(i) + h * model->getAcceleration(i));
            }
        }
    }
}

void SPH::TimeStepPF::matrixFreeRHS(const VectorXr & x, VectorXr & result)
{
    Simulation *sim = Simulation::getCurrent();
    const Real h = TimeManager::getCurrent()->getTimeStepSize();;
    const unsigned int nFluids = sim->numberOfFluidModels();
    const unsigned int nBoundaries = sim->numberOfBoundaryModels();

    // helper functions
    // constraint value
    const auto calculateC = [&](const unsigned int fluidModelIndex, const unsigned int i, std::vector<Vector3r> & p) -> Real
    {
        const FluidModel * model = sim->getFluidModel(fluidModelIndex);
        // Compute current density for particle i
        Real density = model->getVolume(i) * sim->W_zero();
        const Vector3r &xi = p[0];
        unsigned int counter = 1;

        for (unsigned int pid = 0; pid < nFluids; pid++)
        {
            const FluidModel *fm_neighbor = sim->getFluidModelFromPointSet(pid);
            for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
            {
                const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
                const Vector3r & xj = p[counter++];
                density += fm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
            }
        }
        // influence of boundary on density
        if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
        {
            forall_boundary_neighbors(
                // Boundary: Akinci2012
                density += bm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
            );
        }
        else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
        {
            forall_density_maps(
                density += rho;
            );
        }
        else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
        {
            forall_volume_maps(
                density += Vj * sim->W(xi - xj);
            );
        }
        // constraint value = density / density0 - 1
        const auto C = density - 1;
        // pressure clamping
        return (C < 0) ? 0 : C;
    };
    // constraint gradient
    const auto calculateNablaC = [&](const unsigned int fluidModelIndex, const unsigned int i, std::vector<Vector3r> & p) -> std::vector<Vector3r>
    {
        std::vector<Vector3r> nablaC(p.size());
        nablaC[0].setZero();
        const Vector3r &xi = p[0];

        unsigned int counter = 1;
        for (unsigned int pid = 0; pid < nFluids; pid++)
        {
            const FluidModel *fm_neighbor = sim->getFluidModelFromPointSet(pid);
            for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
            {
                const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
                const Vector3r & xj = p[counter];
                nablaC[counter] = -fm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
                nablaC[0] -= nablaC[counter];
                counter++;
            }
        }
        if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
        {
            // influence of boundary on gradient
            forall_boundary_neighbors(
                // Boundary: Akinci2012
                nablaC[0] += bm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj);
            )
        }
        else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
        {
            forall_density_maps(
                nablaC[0] -= gradRho;
            );
        }
        else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
        {
            forall_volume_maps(
                nablaC[0] += Vj * sim->gradW(xi - xj);
            );
        }
        return nablaC;
    };

    AtomicRealVec accumulator(3 * m_numActiveParticlesTotal);

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const unsigned int numParticles = model->numActiveParticles();
        const unsigned int offset = m_simulationData.getParticleOffset(fluidModelIndex);
        const Real density0 = model->getDensity0();

        // influence of pressure
        #pragma omp parallel default(shared)
        {
            // local step for fluid constraints
            #pragma omp for schedule(static)
            for (int i = 0; i < (int)numParticles; i++)
            {
                // total number of neighbors in all point sets (fluids and boundaries)
                unsigned int numParticlesInConstraint = 1;
                for (unsigned int pid = 0; pid < nFluids; pid++)
                {
                    numParticlesInConstraint += sim->numberOfNeighbors(fluidModelIndex, pid, i);
                }
                // particle positions in current constraint, will be projected
                std::vector<Vector3r> p;
                p.reserve(numParticlesInConstraint);
                // the i'th particle itself
                p.emplace_back(x.Vec3Block(offset + i));
                // fluid neighbors
                for (unsigned int pid = 0; pid < nFluids; pid++)
                {
                    const unsigned int neighborOffset = m_simulationData.getParticleOffset(pid);
                    for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
                    {
                        const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
                        p.emplace_back(x.Vec3Block(neighborOffset + neighborIndex));
                    }
                }

                // constraint projection
                const auto C_goal = Real(1e-14);
                const auto max_steps = 100u;
                auto it = 0u;
                auto C = calculateC(fluidModelIndex, i, p);
                while ((std::abs(C) > C_goal) && it++ < max_steps)
                {
                    const auto g = calculateNablaC(fluidModelIndex, i, p);
                    const Real dg = [&g]() { Real s = 0;  for (const auto & x : g) { s += x.squaredNorm(); } return s; }();
                    if (dg == 0) break; // found a minimum
                    const Real cdg = -C / (dg + 1e-6f); // add regularization factor

                    // move fluid particles along constraint gradient
                    p[0] += (cdg * m_simulationData.getNumFluidNeighbors(fluidModelIndex, i)) * g[0];

                    // fluid particles are projected
                    unsigned int counter = 1;
                    for (unsigned int pid = 0; pid < nFluids; pid++)
                    {
                        for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
                        {
                            const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
                            const unsigned int nfn = m_simulationData.getNumFluidNeighbors(pid, neighborIndex);
                            p[counter] += (cdg * nfn) * g[counter];
                            counter++;
                        }
                    }

                    if (it + 1 < max_steps)
                    {
                        if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
                            computeDensityAndGradient(fluidModelIndex, i, p[0]);

                        C = calculateC(fluidModelIndex, i, p);
                    }
                }
                // update RHS
                const unsigned int index_i = offset + i;
                for (auto c = 0u; c < 3; c++)
                    addToAtomicReal(accumulator[3 * index_i + c]._a, density0 * p[0][c]);
                unsigned int counter = 1;
                for (unsigned int pid = 0; pid < nFluids; pid++)
                {
                    const unsigned int neighborOffset = m_simulationData.getParticleOffset(pid);
                    for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
                    {
                        const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
                        const unsigned int index_j = neighborOffset + neighborIndex;
                        for (auto c = 0u; c < 3; c++)
                            addToAtomicReal(accumulator[3 * index_j + c]._a, density0 * p[counter][c]);
                        counter++;
                    }
                }
            }
        }
    }

    const Real system_scale = h * h * m_stiffness;
    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++)
    {
        FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const unsigned int numParticles = model->numActiveParticles();
        const unsigned int offset = m_simulationData.getParticleOffset(fluidModelIndex);
        const Real density_scaled_system = system_scale / model->getDensity0();

        // influence of momentum
        #pragma omp parallel default(shared)
        {
            #pragma omp for schedule(static)
            for (int i = 0; i < (int)numParticles; i++)
            {
                const auto  m = model->getMass(i);
                const auto& s = m_simulationData.getS(fluidModelIndex, i);
                for (auto c = 0u; c < 3; c++)
                {
                    const auto id = 3 * (offset + i) + c;
                    result[id] = density_scaled_system * accumulator[id]._a + m * s[c];
                }
            }
        }
    }
}

void SPH::TimeStepPF::matrixVecProd(const Real* vec, Real *result, void *userData)
{
    TimeStepPF *timeStepPF = static_cast<TimeStepPF*>(userData);

    Simulation *sim = Simulation::getCurrent();
    const Real  h = TimeManager::getCurrent()->getTimeStepSize();
    const unsigned int nFluids = sim->numberOfFluidModels();

    const Real system_scale = h * h * timeStepPF->m_stiffness;

    for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++)
    {
        const FluidModel *model = sim->getFluidModel(fluidModelIndex);
        const Real density0 = model->getDensity0();
        const Real density_scaled_system = system_scale / density0;
        const unsigned int numParticles = model->numActiveParticles();
        const unsigned int offset = timeStepPF->m_simulationData.getParticleOffset(fluidModelIndex);

        #pragma omp parallel for schedule(static)
        for (int i = 0; i < (int) numParticles; i++)
        {
            const unsigned int index_i = i + offset;

            Eigen::Map<Vector3r> ri(result + 3 * index_i, 3, 1);
            Eigen::Map<const Vector3r> xi(vec + 3 * index_i, 3, 1);

            // the particle itself
            Real accumulator = density0;
            // fluid neighbors
            for (unsigned int pid = 0; pid < nFluids; pid++)
            {
                const FluidModel *fm_neighbor = sim->getFluidModelFromPointSet(pid);
                accumulator += fm_neighbor->getDensity0() * sim->numberOfNeighbors(fluidModelIndex, pid, i);
            }
            // boundary neighbors are not part of the linear system,
            // they just influence the right hand side in the projection

            ri = (density_scaled_system * accumulator + model->getMass(i)) * xi;
        }
    }
}

void TimeStepPF::performNeighborhoodSearchSort()
{
    m_simulationData.performNeighborhoodSearchSort();
}

void TimeStepPF::emittedParticles(FluidModel *model, const unsigned int startIndex)
{
    m_simulationData.emittedParticles(model, startIndex);
}

void TimeStepPF::resize()
{
    m_simulationData.init();
}