.. _program_listing_file_SPlisHSPlasH_Simulation.cpp:

Program Listing for File Simulation.cpp
=======================================

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

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

.. code-block:: cpp

   #include "Simulation.h"
   #include "TimeManager.h"
   #include "Utilities/Timing.h"
   #include "TimeStep.h"
   #include "EmitterSystem.h"
   #include "SPlisHSPlasH/WCSPH/TimeStepWCSPH.h"
   #include "SPlisHSPlasH/PCISPH/TimeStepPCISPH.h"
   #include "SPlisHSPlasH/PBF/TimeStepPBF.h"
   #include "SPlisHSPlasH/IISPH/TimeStepIISPH.h"
   #include "SPlisHSPlasH/DFSPH/TimeStepDFSPH.h"
   #include "SPlisHSPlasH/PF/TimeStepPF.h"
   #include "SPlisHSPlasH/ICSPH/TimeStepICSPH.h"
   #include "BoundaryModel_Akinci2012.h"
   #include "BoundaryModel_Bender2019.h"
   #include "BoundaryModel_Koschier2017.h"
   
   
   using namespace SPH;
   using namespace std;
   using namespace GenParam;
   
   Simulation* Simulation::current = nullptr;
   int Simulation::SIM_2D = -1;
   int Simulation::PARTICLE_RADIUS = -1;
   int Simulation::GRAVITATION = -1;
   int Simulation::CFL_METHOD = -1;
   int Simulation::CFL_FACTOR = -1;
   int Simulation::CFL_MIN_TIMESTEPSIZE = -1;
   int Simulation::CFL_MAX_TIMESTEPSIZE = -1;
   int Simulation::ENABLE_Z_SORT = -1;
   int Simulation::STEPS_PER_Z_SORT = -1;
   int Simulation::KERNEL_METHOD = -1;
   int Simulation::GRAD_KERNEL_METHOD = -1;
   int Simulation::ENUM_KERNEL_CUBIC = -1;
   int Simulation::ENUM_KERNEL_WENDLANDQUINTICC2 = -1;
   int Simulation::ENUM_KERNEL_POLY6 = -1;
   int Simulation::ENUM_KERNEL_SPIKY = -1;
   int Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC = -1;
   int Simulation::ENUM_KERNEL_CUBIC_2D = -1;
   int Simulation::ENUM_KERNEL_WENDLANDQUINTICC2_2D = -1;
   int Simulation::ENUM_GRADKERNEL_CUBIC = -1;
   int Simulation::ENUM_GRADKERNEL_WENDLANDQUINTICC2 = -1;
   int Simulation::ENUM_GRADKERNEL_POLY6 = -1;
   int Simulation::ENUM_GRADKERNEL_SPIKY = -1;
   int Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC = -1;
   int Simulation::ENUM_GRADKERNEL_CUBIC_2D = -1;
   int Simulation::ENUM_GRADKERNEL_WENDLANDQUINTICC2_2D = -1;
   int Simulation::SIMULATION_METHOD = -1;
   int Simulation::ENUM_CFL_NONE = -1;
   int Simulation::ENUM_CFL_STANDARD = -1;
   int Simulation::ENUM_CFL_ITER = -1;
   int Simulation::ENUM_SIMULATION_WCSPH = -1;
   int Simulation::ENUM_SIMULATION_PCISPH = -1;
   int Simulation::ENUM_SIMULATION_PBF = -1;
   int Simulation::ENUM_SIMULATION_IISPH = -1;
   int Simulation::ENUM_SIMULATION_DFSPH = -1;
   int Simulation::ENUM_SIMULATION_PF = -1;
   int Simulation::ENUM_SIMULATION_ICSPH = -1;
   int Simulation::BOUNDARY_HANDLING_METHOD = -1;
   int Simulation::ENUM_AKINCI2012 = -1;
   int Simulation::ENUM_KOSCHIER2017 = -1;
   int Simulation::ENUM_BENDER2019 = -1;
   
   
   Simulation::Simulation () 
   {
       m_cflMethod = 1;
       m_cflFactor = 0.5;
       m_cflMinTimeStepSize = static_cast<Real>(0.0001);
       m_cflMaxTimeStepSize = static_cast<Real>(0.005);
       m_gravitation = Vector3r(0.0, static_cast<Real>(-9.81), 0.0);
   
       m_kernelMethod = -1;
       m_gradKernelMethod = -1;
   
       m_neighborhoodSearch = nullptr;
       m_timeStep = nullptr;
       m_simulationMethod = SimulationMethods::NumSimulationMethods;
       m_simulationMethodChanged = NULL;
   
       m_simulationIsInitialized = false;
       m_sim2D = false;
       m_enableZSort = true;
       m_stepsPerZSort = 500;
       m_counter = 0;
   
       m_animationFieldSystem = new AnimationFieldSystem();
       m_boundaryHandlingMethod = static_cast<int>(BoundaryHandlingMethods::Bender2019);
   }
   
   Simulation::~Simulation () 
   {
   #ifdef USE_DEBUG_TOOLS
       delete m_debugTools;
   #endif
       delete m_animationFieldSystem;
       delete m_timeStep;
       delete m_neighborhoodSearch;
       delete TimeManager::getCurrent();
   
       for (unsigned int i = 0; i < m_fluidModels.size(); i++)
           delete m_fluidModels[i];
       m_fluidModels.clear();
   
       for (unsigned int i = 0; i < m_boundaryModels.size(); i++)
           delete m_boundaryModels[i];
       m_boundaryModels.clear();
   
       current = nullptr;
   }
   
   Simulation* Simulation::getCurrent ()
   {
       if (current == nullptr)
       {
           current = new Simulation ();
       }
       return current;
   }
   
   void Simulation::setCurrent (Simulation* tm)
   {
       current = tm;
   }
   
   bool Simulation::hasCurrent()
   {
       return (current != nullptr);
   }
   
   void Simulation::init(const Real particleRadius, const bool sim2D)
   {
       m_sim2D = sim2D;
       initParameters();
   
       registerNonpressureForces();
   
       // init kernel
       setParticleRadius(particleRadius);
   
       setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC);
       setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
   
       // Initialize neighborhood search
       if (m_neighborhoodSearch == NULL)
   #ifdef GPU_NEIGHBORHOOD_SEARCH
           m_neighborhoodSearch = new NeighborhoodSearch(m_supportRadius);
   #else
           m_neighborhoodSearch = new NeighborhoodSearch(m_supportRadius, false);
   #endif
       m_neighborhoodSearch->set_radius(m_supportRadius);
   }
   
   void Simulation::deferredInit()
   {
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
       {
           FluidModel* fm = getFluidModel(i);
           fm->deferredInit();
       }
   }
   
   void Simulation::initParameters()
   {
       ParameterObject::initParameters();
   
       SIM_2D = createBoolParameter("sim2D", "2D Simulation", &m_sim2D);
       setGroup(SIM_2D, "Simulation|Simulation");
       setDescription(SIM_2D, "2D/3D simulation.");
       getParameter(SIM_2D)->setReadOnly(true);
   
       ENABLE_Z_SORT = createBoolParameter("enableZSort", "Enable z-sort", &m_enableZSort);
       setGroup(ENABLE_Z_SORT, "Simulation|Simulation");
       setDescription(ENABLE_Z_SORT, "Enable z-sort to improve cache hits.");
   
       STEPS_PER_Z_SORT = createNumericParameter("stepsPerZSort", "Simulation steps per z-sort", &m_stepsPerZSort);
       setGroup(STEPS_PER_Z_SORT, "Simulation|Simulation");
       setDescription(STEPS_PER_Z_SORT, "Number of simulation steps which are performed before a z-sort is applied.");
       static_cast<UnsignedIntParameter*>(getParameter(STEPS_PER_Z_SORT))->setMinValue(1u);
   
       ParameterBase::GetFunc<Real> getRadiusFct = std::bind(&Simulation::getParticleRadius, this);
       ParameterBase::SetFunc<Real> setRadiusFct = std::bind(&Simulation::setParticleRadius, this, std::placeholders::_1);
       PARTICLE_RADIUS = createNumericParameter("particleRadius", "Particle radius", getRadiusFct, setRadiusFct);
       setGroup(PARTICLE_RADIUS, "Simulation|Simulation");
       setDescription(PARTICLE_RADIUS, "Radius of the fluid particles.");
       getParameter(PARTICLE_RADIUS)->setReadOnly(true);
   
       GRAVITATION = createVectorParameter("gravitation", "Gravitation", 3u, m_gravitation.data());
       setGroup(GRAVITATION, "Simulation|Simulation");
       setDescription(GRAVITATION, "Vector to define the gravitational acceleration.");
   
       CFL_METHOD = createEnumParameter("cflMethod", "CFL - method", &m_cflMethod);
       setGroup(CFL_METHOD, "Simulation|CFL");
       setDescription(CFL_METHOD, "CFL method used for adaptive time stepping.");
       EnumParameter *enumParam = static_cast<EnumParameter*>(getParameter(CFL_METHOD));
       enumParam->addEnumValue("None", ENUM_CFL_NONE);
       enumParam->addEnumValue("CFL", ENUM_CFL_STANDARD);
       enumParam->addEnumValue("CFL - iterations", ENUM_CFL_ITER);
   
       CFL_FACTOR = createNumericParameter("cflFactor", "CFL - factor", &m_cflFactor);
       setGroup(CFL_FACTOR, "Simulation|CFL");
       setDescription(CFL_FACTOR, "Factor to scale the CFL time step size.");
       static_cast<RealParameter*>(getParameter(CFL_FACTOR))->setMinValue(static_cast<Real>(1e-6));
   
       CFL_MIN_TIMESTEPSIZE = createNumericParameter("cflMinTimeStepSize", "CFL - min. time step size", &m_cflMinTimeStepSize);
       setGroup(CFL_MIN_TIMESTEPSIZE, "Simulation|CFL");
       setDescription(CFL_MIN_TIMESTEPSIZE, "Min. time step size.");
       static_cast<RealParameter*>(getParameter(CFL_MIN_TIMESTEPSIZE))->setMinValue(static_cast<Real>(1e-9));
   
       CFL_MAX_TIMESTEPSIZE = createNumericParameter("cflMaxTimeStepSize", "CFL - max. time step size", &m_cflMaxTimeStepSize);
       setGroup(CFL_MAX_TIMESTEPSIZE, "Simulation|CFL");
       setDescription(CFL_MAX_TIMESTEPSIZE, "Max. time step size.");
       static_cast<RealParameter*>(getParameter(CFL_MAX_TIMESTEPSIZE))->setMinValue(static_cast<Real>(1e-6));
   
       ParameterBase::GetFunc<int> getKernelFct = std::bind(&Simulation::getKernel, this);
       ParameterBase::SetFunc<int> setKernelFct = std::bind(&Simulation::setKernel, this, std::placeholders::_1);
       KERNEL_METHOD = createEnumParameter("kernel", "Kernel", getKernelFct, setKernelFct);
       setGroup(KERNEL_METHOD, "Simulation|Kernel");
       setDescription(KERNEL_METHOD, "Kernel function used in the SPH model.");
       enumParam = static_cast<EnumParameter*>(getParameter(KERNEL_METHOD));
       if (!m_sim2D)
       {
           enumParam->addEnumValue("Cubic spline", ENUM_KERNEL_CUBIC);
           enumParam->addEnumValue("Wendland quintic C2", ENUM_KERNEL_WENDLANDQUINTICC2);
           enumParam->addEnumValue("Poly6", ENUM_KERNEL_POLY6);
           enumParam->addEnumValue("Spiky", ENUM_KERNEL_SPIKY);
           enumParam->addEnumValue("Precomputed cubic spline", ENUM_KERNEL_PRECOMPUTED_CUBIC);
       }
       else
       {
           enumParam->addEnumValue("Cubic spline 2D", ENUM_KERNEL_CUBIC_2D);
           enumParam->addEnumValue("Wendland quintic C2 2D", ENUM_KERNEL_WENDLANDQUINTICC2_2D);
       }
   
       ParameterBase::GetFunc<int> getGradKernelFct = std::bind(&Simulation::getGradKernel, this);
       ParameterBase::SetFunc<int> setGradKernelFct = std::bind(&Simulation::setGradKernel, this, std::placeholders::_1);
       GRAD_KERNEL_METHOD = createEnumParameter("gradKernel", "Gradient of kernel", getGradKernelFct, setGradKernelFct);
       setGroup(GRAD_KERNEL_METHOD, "Simulation|Kernel");
       setDescription(GRAD_KERNEL_METHOD, "Gradient of the kernel function used in the SPH model.");
       enumParam = static_cast<EnumParameter*>(getParameter(GRAD_KERNEL_METHOD));
       if (!m_sim2D)
       {
           enumParam->addEnumValue("Cubic spline", ENUM_GRADKERNEL_CUBIC);
           enumParam->addEnumValue("Wendland quintic C2", ENUM_GRADKERNEL_WENDLANDQUINTICC2);
           enumParam->addEnumValue("Poly6", ENUM_GRADKERNEL_POLY6);
           enumParam->addEnumValue("Spiky", ENUM_GRADKERNEL_SPIKY);
           enumParam->addEnumValue("Precomputed cubic spline", ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
       }
       else
       {
           enumParam->addEnumValue("Cubic spline 2D", ENUM_GRADKERNEL_CUBIC_2D);
           enumParam->addEnumValue("Wendland quintic C2 2D", ENUM_GRADKERNEL_WENDLANDQUINTICC2_2D);
       }
   
       ParameterBase::GetFunc<int> getSimulationFct = std::bind(&Simulation::getSimulationMethod, this);
       ParameterBase::SetFunc<int> setSimulationFct = std::bind(&Simulation::setSimulationMethod, this, std::placeholders::_1);
       SIMULATION_METHOD = createEnumParameter("simulationMethod", "Simulation method", getSimulationFct, setSimulationFct);
       setGroup(SIMULATION_METHOD, "Simulation|Simulation");
       setDescription(SIMULATION_METHOD, "Simulation method.");
       enumParam = static_cast<EnumParameter*>(getParameter(SIMULATION_METHOD));
       enumParam->addEnumValue(TimeStepWCSPH::METHOD_NAME, ENUM_SIMULATION_WCSPH);
       enumParam->addEnumValue(TimeStepPCISPH::METHOD_NAME, ENUM_SIMULATION_PCISPH);
       enumParam->addEnumValue(TimeStepPBF::METHOD_NAME, ENUM_SIMULATION_PBF);
       enumParam->addEnumValue(TimeStepIISPH::METHOD_NAME, ENUM_SIMULATION_IISPH);
       enumParam->addEnumValue(TimeStepDFSPH::METHOD_NAME, ENUM_SIMULATION_DFSPH);
       enumParam->addEnumValue(TimeStepPF::METHOD_NAME, ENUM_SIMULATION_PF);
       enumParam->addEnumValue(TimeStepICSPH::METHOD_NAME, ENUM_SIMULATION_ICSPH);
   
       BOUNDARY_HANDLING_METHOD = createEnumParameter("boundaryHandlingMethod", "Boundary handling method", &m_boundaryHandlingMethod);
       setGroup(BOUNDARY_HANDLING_METHOD, "Simulation|Simulation");
       setDescription(BOUNDARY_HANDLING_METHOD, "Boundary handling method.");
       enumParam = static_cast<EnumParameter*>(getParameter(BOUNDARY_HANDLING_METHOD));
       enumParam->addEnumValue("Akinci et al. 2012", ENUM_AKINCI2012);
       enumParam->addEnumValue("Koschier and Bender 2017", ENUM_KOSCHIER2017);
       enumParam->addEnumValue("Bender et al. 2019", ENUM_BENDER2019);
       enumParam->setReadOnly(true);
   }
   
   
   void Simulation::setParticleRadius(Real val)
   {
       m_particleRadius = val;
       m_supportRadius = static_cast<Real>(4.0) * m_particleRadius;
       initKernels();
   }
   
   void Simulation::initKernels()
   {
       // init kernel
       Poly6Kernel::setRadius(m_supportRadius);
       SpikyKernel::setRadius(m_supportRadius);
       CubicKernel::setRadius(m_supportRadius);
       WendlandQuinticC2Kernel::setRadius(m_supportRadius);
       PrecomputedCubicKernel::setRadius(m_supportRadius);
       CohesionKernel::setRadius(m_supportRadius);
       AdhesionKernel::setRadius(m_supportRadius);
       CubicKernel2D::setRadius(m_supportRadius);
       WendlandQuinticC2Kernel2D::setRadius(m_supportRadius);
   #ifdef USE_AVX
       CubicKernel_AVX::setRadius(m_supportRadius, m_sim2D);
       //  Poly6Kernel_AVX::setRadius(m_supportRadius);
       //  SpikyKernel_AVX::setRadius(m_supportRadius);
   #endif
   }
   
   void Simulation::setGradKernel(int val)
   {
       m_gradKernelMethod = val;
   
       if (!m_sim2D)
       {
           if ((m_gradKernelMethod < 0) || (m_gradKernelMethod > 4))
               m_gradKernelMethod = 0;
   
           if (m_gradKernelMethod == 0)
               m_gradKernelFct = CubicKernel::gradW;
           else if (m_gradKernelMethod == 1)
               m_gradKernelFct = WendlandQuinticC2Kernel::gradW;
           else if (m_gradKernelMethod == 2)
               m_gradKernelFct = Poly6Kernel::gradW;
           else if (m_gradKernelMethod == 3)
               m_gradKernelFct = SpikyKernel::gradW;
           else if (m_gradKernelMethod == 4)
               m_gradKernelFct = Simulation::PrecomputedCubicKernel::gradW;
       }
       else
       {
           if ((m_gradKernelMethod < 0) || (m_gradKernelMethod > 1))
               m_gradKernelMethod = 0;
   
           if (m_gradKernelMethod == 0)
               m_gradKernelFct = CubicKernel2D::gradW;
           else if (m_gradKernelMethod == 1)
               m_gradKernelFct = WendlandQuinticC2Kernel2D::gradW;
       }
   }
   
   void Simulation::setKernel(int val)
   {
       if (val == m_kernelMethod)
           return;
   
       m_kernelMethod = val;
       if (!m_sim2D)
       {
           if ((m_kernelMethod < 0) || (m_kernelMethod > 4))
               m_kernelMethod = 0;
   
           if (m_kernelMethod == 0)
           {
               m_W_zero = CubicKernel::W_zero();
               m_kernelFct = CubicKernel::W;
           }
           else if (m_kernelMethod == 1)
           {
               m_W_zero = WendlandQuinticC2Kernel::W_zero();
               m_kernelFct = WendlandQuinticC2Kernel::W;
           }
           else if (m_kernelMethod == 2)
           {
               m_W_zero = Poly6Kernel::W_zero();
               m_kernelFct = Poly6Kernel::W;
           }
           else if (m_kernelMethod == 3)
           {
               m_W_zero = SpikyKernel::W_zero();
               m_kernelFct = SpikyKernel::W;
           }
           else if (m_kernelMethod == 4)
           {
               m_W_zero = Simulation::PrecomputedCubicKernel::W_zero();
               m_kernelFct = Simulation::PrecomputedCubicKernel::W;
           }
       }
       else
       {
           if ((m_kernelMethod < 0) || (m_kernelMethod > 1))
               m_kernelMethod = 0;
   
           if (m_kernelMethod == 0)
           {
               m_W_zero = CubicKernel2D::W_zero();
               m_kernelFct = CubicKernel2D::W;
           }
           else if (m_kernelMethod == 1)
           {
               m_W_zero = WendlandQuinticC2Kernel2D::W_zero();
               m_kernelFct = WendlandQuinticC2Kernel2D::W;
           }
       }
       if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
           updateBoundaryVolume();
   }
   
   void Simulation::updateTimeStepSize()
   {
       if (m_cflMethod == 1)
           updateTimeStepSizeCFL();
       else if (m_cflMethod == 2)
       {
           Real h = TimeManager::getCurrent()->getTimeStepSize();
           updateTimeStepSizeCFL();
           const unsigned int iterations = m_timeStep->getNumIterations();
           if (iterations == 0)        // check if the method requried iterations
               return;
           if (iterations > 10)
               h *= static_cast<Real>(0.9);
           else if (iterations < 5)
               h *= static_cast<Real>(1.1);
           h = min(h, TimeManager::getCurrent()->getTimeStepSize());
           TimeManager::getCurrent()->setTimeStepSize(h);
       }
   }
   
   void Simulation::updateTimeStepSizeCFL()
   {
       const Real radius = m_particleRadius;
       Real h = TimeManager::getCurrent()->getTimeStepSize();
       Simulation *sim = Simulation::getCurrent();
       const unsigned int nBoundaries = sim->numberOfBoundaryModels();
   
       // Approximate max. position change due to current velocities
       Real maxVel = 0.0;
       const Real diameter = static_cast<Real>(2.0)*radius;
   
       // fluid particles
       for (unsigned int fluidModelIndex = 0; fluidModelIndex < numberOfFluidModels(); fluidModelIndex++)
       {
           FluidModel *fm = getFluidModel(fluidModelIndex);
           const unsigned int numParticles = fm->numActiveParticles();
           for (unsigned int i = 0; i < numParticles; i++)
           {
               const Vector3r &vel = fm->getVelocity(i);
               const Vector3r &accel = fm->getAcceleration(i);
               const Real velMag = (vel + accel*h).squaredNorm();
               if (velMag > maxVel)
                   maxVel = velMag;
           }
       }
   
       // boundary particles
       if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
       {
           for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
           {
               BoundaryModel_Akinci2012 *bm = static_cast<BoundaryModel_Akinci2012*>(getBoundaryModel(i));
               if (bm->getRigidBodyObject()->isDynamic() || bm->getRigidBodyObject()->isAnimated())
               {
                   for (unsigned int j = 0; j < bm->numberOfParticles(); j++)
                   {
                       const Vector3r &vel = bm->getVelocity(j);
                       const Real velMag = vel.squaredNorm();
                       if (velMag > maxVel)
                           maxVel = velMag;
                   }
               }
           }
       }
       else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
       {
           for (unsigned int boundaryModelIndex = 0; boundaryModelIndex < numberOfBoundaryModels(); boundaryModelIndex++)
           {
               BoundaryModel_Koschier2017 *bm = static_cast<BoundaryModel_Koschier2017*>(getBoundaryModel(boundaryModelIndex));
               if (bm->getRigidBodyObject()->isDynamic() || bm->getRigidBodyObject()->isAnimated())
               {
                   maxVel = std::max(maxVel, bm->getMaxVel());
               }
           }
       }
       else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
       {
           for (unsigned int boundaryModelIndex = 0; boundaryModelIndex < numberOfBoundaryModels(); boundaryModelIndex++)
           {
               BoundaryModel_Bender2019 *bm = static_cast<BoundaryModel_Bender2019*>(getBoundaryModel(boundaryModelIndex));
               if (bm->getRigidBodyObject()->isDynamic() || bm->getRigidBodyObject()->isAnimated())
               {
                   maxVel = std::max(maxVel, bm->getMaxVel());
               }
           }
       }
   
       // avoid division by zero
       if (maxVel < static_cast<Real>(1.0e-9))
           maxVel = static_cast<Real>(1.0e-9);
   
       // Approximate max. time step size      
       h = m_cflFactor * static_cast<Real>(0.4) * (diameter / (sqrt(maxVel)));
   
       h = min(h, m_cflMaxTimeStepSize);
       h = max(h, m_cflMinTimeStepSize);
   
       TimeManager::getCurrent()->setTimeStepSize(h);
   }
   
   void Simulation::computeNonPressureForces()
   {
       START_TIMING("computeNonPressureForces")
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
       {
           FluidModel *fm = getFluidModel(i);
           fm->computeVorticity();
           fm->computeSurfaceTension();
           fm->computeViscosity();
           fm->computeDragForce();
           fm->computeElasticity();
           fm->computeXSPH();
       }
       STOP_TIMING_AVG
   }
   
   void Simulation::reset()
   {
       // reset fluid models
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
           getFluidModel(i)->reset();
   
       // reset boundary models
       for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
           getBoundaryModel(i)->reset();
   
       if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
           updateBoundaryVolume();
   
       performNeighborhoodSearchSort();
       if (m_timeStep)
           m_timeStep->reset();
   
       m_animationFieldSystem->reset();
   
       m_counter = 0;
   
       TimeManager::getCurrent()->setTime(0.0);
   }
   
   void Simulation::setSimulationMethod(const int val)
   {
       SimulationMethods method = static_cast<SimulationMethods>(val);
       if ((method < SimulationMethods::WCSPH) || (method >= SimulationMethods::NumSimulationMethods))
           method = SimulationMethods::DFSPH;
   
       if (method == m_simulationMethod)
           return;
   
       delete m_timeStep;
       m_timeStep = nullptr;
   
       m_simulationMethod = method;
   
       if (method == SimulationMethods::WCSPH)
       {
           m_timeStep = new TimeStepWCSPH();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_CUBIC);
       }
       else if (method == SimulationMethods::PCISPH)
       {
           m_timeStep = new TimeStepPCISPH();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_CUBIC);
       }
       else if (method == SimulationMethods::PBF)
       {
           m_timeStep = new TimeStepPBF();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_POLY6);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_SPIKY);
       }
       else if (method == SimulationMethods::IISPH)
       {
           m_timeStep = new TimeStepIISPH();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_CUBIC);
       }
       else if (method == SimulationMethods::DFSPH)
       {
           m_timeStep = new TimeStepDFSPH();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
       }
       else if (method == SimulationMethods::PF)
       {
           m_timeStep = new TimeStepPF();
           m_timeStep->init();
   
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
       }
       else if (method == SimulationMethods::ICSPH)
       {
           m_timeStep = new TimeStepICSPH();
           m_timeStep->init();
           setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_CUBIC);
           setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_CUBIC);
       }
   
       if (m_simulationMethodChanged != nullptr)
           m_simulationMethodChanged();
   }
   
   
   
   void Simulation::performNeighborhoodSearch()
   {
       if (zSortEnabled())
       {
           if (m_counter % m_stepsPerZSort == 0)
           {
               performNeighborhoodSearchSort();
           }
           m_counter++;
       }
       START_TIMING("neighborhood_search");
       m_neighborhoodSearch->find_neighbors();
       STOP_TIMING_AVG;
   }
   
   void Simulation::performNeighborhoodSearchSort()
   {
       if (!zSortEnabled())
           return;
   
       m_neighborhoodSearch->z_sort();
   
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
       {
           FluidModel *fm = getFluidModel(i);
           fm->performNeighborhoodSearchSort();
       }
       for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
       {
           BoundaryModel *bm = getBoundaryModel(i);
           bm->performNeighborhoodSearchSort();
       }
       getTimeStep()->performNeighborhoodSearchSort();
   #ifdef USE_DEBUG_TOOLS
       m_debugTools->performNeighborhoodSearchSort();
   #endif
   }
   
   void Simulation::setSimulationMethodChangedCallback(std::function<void()> const& callBackFct)
   {
       m_simulationMethodChanged = callBackFct;
   }
   
   void Simulation::setSimulationInitialized(int val) 
   { 
       m_simulationIsInitialized = val; 
       if (m_simulationIsInitialized)
           deferredInit();
   }
   
   void Simulation::emittedParticles(FluidModel *model, const unsigned int startIndex)
   {
       model->emittedParticles(startIndex);
       m_timeStep->emittedParticles(model, startIndex);
   #ifdef USE_DEBUG_TOOLS
       m_debugTools->emittedParticles(model, startIndex);
   #endif
   }
   
   void Simulation::emitParticles()
   {
       START_TIMING("emitParticles");
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
       {
           FluidModel *fm = getFluidModel(i);
           fm->getEmitterSystem()->step();
       }
       STOP_TIMING_AVG
   }
   
   void Simulation::animateParticles()
   {
       START_TIMING("animateParticles");
       m_animationFieldSystem->step();
       STOP_TIMING_AVG
   }
   
   void Simulation::addBoundaryModel(BoundaryModel *bm)
   {
       m_boundaryModels.push_back(bm);
   }
   
   void Simulation::addFluidModel(const std::string &id, const unsigned int nFluidParticles, Vector3r* fluidParticles, Vector3r* fluidVelocities, unsigned int* fluidObjectIds, const unsigned int nMaxEmitterParticles)
   {
       FluidModel *fm = new FluidModel();
       fm->initModel(id, nFluidParticles, fluidParticles, fluidVelocities, fluidObjectIds, nMaxEmitterParticles);
       m_fluidModels.push_back(fm);
   }
   
   
   void Simulation::updateBoundaryVolume()
   {
       if (m_neighborhoodSearch == nullptr)
           return;
   
       Simulation *sim = Simulation::getCurrent();
       const unsigned int nFluids = sim->numberOfFluidModels();
   
       // Compute value psi for boundary particles (boundary handling)
       // (see Akinci et al. "Versatile rigid - fluid coupling for incompressible SPH", Siggraph 2012
   
       // Search boundary neighborhood
   
       // Activate only static boundaries
       LOG_INFO << "Initialize boundary volume";
       m_neighborhoodSearch->set_active(false);
       for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
       {
           if (!getBoundaryModel(i)->getRigidBodyObject()->isDynamic() && !getBoundaryModel(i)->getRigidBodyObject()->isAnimated())
               m_neighborhoodSearch->set_active(i + nFluids, true, true);
       }
   
       //performNeighborhoodSearchSort();
       m_neighborhoodSearch->find_neighbors();
   
       // Boundary objects
       for (unsigned int body = 0; body < numberOfBoundaryModels(); body++)
       {
           if (!getBoundaryModel(body)->getRigidBodyObject()->isDynamic() && !getBoundaryModel(body)->getRigidBodyObject()->isAnimated())
               static_cast<BoundaryModel_Akinci2012*>(getBoundaryModel(body))->computeBoundaryVolume();
       }
   
       // Compute boundary psi for all dynamic bodies
       for (unsigned int body = 0; body < numberOfBoundaryModels(); body++)
       {
           // Deactivate all
           m_neighborhoodSearch->set_active(false);
   
           // Only activate next dynamic body
           if (getBoundaryModel(body)->getRigidBodyObject()->isDynamic() || getBoundaryModel(body)->getRigidBodyObject()->isAnimated())
           {
               m_neighborhoodSearch->set_active(body + nFluids, true, true);
               m_neighborhoodSearch->find_neighbors();
               static_cast<BoundaryModel_Akinci2012*>(getBoundaryModel(body))->computeBoundaryVolume();
           }
       }
   
       // Activate only fluids 
       m_neighborhoodSearch->set_active(false);
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
       {
           for (unsigned int j = 0; j < numberOfFluidModels(); j++)
               m_neighborhoodSearch->set_active(i, j, true);
           for (unsigned int j = numberOfFluidModels(); j < m_neighborhoodSearch->point_sets().size(); j++)
               m_neighborhoodSearch->set_active(i, j, true);
       }
   }
   
   void SPH::Simulation::saveState(BinaryFileWriter &binWriter)
   {
       binWriter.write(m_W_zero);
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
           getFluidModel(i)->saveState(binWriter);
       for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
           getBoundaryModel(i)->saveState(binWriter);
       m_timeStep->saveState(binWriter);
   }
   
   void SPH::Simulation::loadState(BinaryFileReader &binReader)
   {
       binReader.read(m_W_zero);
       for (unsigned int i = 0; i < numberOfFluidModels(); i++)
           getFluidModel(i)->loadState(binReader);
       for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
           getBoundaryModel(i)->loadState(binReader);
       m_timeStep->loadState(binReader);
   }