src/Solver.cpp

Go to the documentation of this file.

#include "Solver.h"
#include "Configuration.h"
#include "Integration.h"

#include "boost/foreach.hpp"

#include "ctime"

Solver::Solver()
{
    //create particles and initialise parameters
    Configuration::initialiseFluidSolver
            (
                m_smoothingLength,
                m_particleList,
                m_hoseParticlePrototypeList,
                m_hoseParticleList,
                m_hoseCenter,
                m_hoseVelocity,
                m_drawHoseMarker,
                m_hoseWaitUntilHitBoundary,
                m_hoseWaitUntilHitRBD
           );

    //precompute kernel weight multipliers
    m_weightPoly = (315.0 / (64.0 * ngl::PI * pow(m_smoothingLength, 9)));
    m_weightPolyGradient = (-945.0 / (32.0 * ngl::PI * pow(m_smoothingLength, 9)));
    m_weightPolyLaplacian = (-945.0 / (32.0 * ngl::PI * pow(m_smoothingLength, 9)));
    m_weightPressureGradient = (-45.0 / (ngl::PI * pow(m_smoothingLength, 6)));
    m_weightViscosityLaplacian = (45.0 / (ngl::PI * pow(m_smoothingLength, 6)));

    //initialise parameters
    m_enableAutoFluid = true;
    m_currentHoseableFluid = 0;

    //create neighbouring structure
    m_neighbour = new Neighbour(m_particleList.size(), m_smoothingLength);

    m_iteration = 0;
}

Solver::~Solver()
{
    //clean up
    std::cout << "Fluid Solver Cleanup" << std::endl;

    //delete objects
    delete m_neighbour;
}

void Solver::initialiseDrawing()
{
    ngl::VBOPrimitives* prim = ngl::VBOPrimitives::instance();

    //create vbo for particle sphere
    prim->createVBOSphere("particleSphere", 0.17, 10);
}

ngl::Real Solver::WKernel_Poly(const ngl::Real _separationLength)
{
    ngl::Real weight = m_weightPoly * (pow(pow(m_smoothingLength, 2) - pow(_separationLength, 2), 3));

    return weight;
}

ngl::Vector Solver::WKernel_Poly_Gradient(const ngl::Vector _separationVector)
{
    ngl::Vector weight = 0;
    ngl::Real separationLength = _separationVector.length();

    weight = m_weightPolyGradient * _separationVector * pow(pow(m_smoothingLength, 2) - pow(separationLength, 2), 2);

    return weight;
}

ngl::Real Solver::WKernel_Poly_Laplacian(const ngl::Real _separationLength)
{
    ngl::Real weight = m_weightPolyLaplacian * (pow(m_smoothingLength, 2) - pow(_separationLength, 2)) * ((3.0 * pow(m_smoothingLength, 2)) - (7.0 * pow(_separationLength, 2)));

    return weight;
}

ngl::Vector Solver::WKernel_Pressure_Gradient(const ngl::Vector _separationVector)
{
    ngl::Vector weight = 0;
    ngl::Real separationLength = _separationVector.length();

    weight = m_weightPressureGradient * pow(m_smoothingLength - separationLength, 2) * normalise(_separationVector);

    return weight;
}

ngl::Real Solver::WKernel_Viscosity_Laplacian(const ngl::Real _separationLength)
{
    ngl::Real weight = m_weightViscosityLaplacian * (m_smoothingLength - _separationLength);

    return weight;
}

ngl::Vector Solver::normalise(const ngl::Vector _vector)
{
    ngl::Vector v = _vector;
    if (_vector.length() > 0.0001f) v.normalize();

    return v;
}


int Solver::getFluidIdFromName(const std::string _name)
{
    //assume id to be of the first fluid
    int id = 0;

    //loop through hose prototype list
    for (int i = 0; i < m_hoseParticlePrototypeList.size(); i++)
    {
        //compare names
        if (m_hoseParticlePrototypeList[i].getName() == _name)
        {
            id = i;
            break;
        }
    }

    return id;
}

void Solver::injectParticles()
{
    std::vector<FluidParticle> l;

    //get translate in hose center
    ngl::Vector translatedVector = m_hoseCenter - m_hoseInitialCenter;

    //get prototype particle from hose prototype list
    FluidParticle prototypeParticle = m_hoseParticlePrototypeList[m_currentHoseableFluid];

    //loop through hose particles and add particles to solver particle list
    for (int i = 0; i < m_hoseParticleList.size(); i++)
    {
        //create particle from prototype particle
        FluidParticle p(prototypeParticle);

        //set id of new particle
        p.setId(Configuration::s_nextParticleId++);

        //set position of new particle = hose particle(i) + translatedVector
        p.setPosition(m_hoseParticleList[i].getPosition() + translatedVector);

        //set velocity of new particle
        p.setVelocity(m_hoseVelocity);

        //set wait until first hit properties
        p.setWaitUntilFirstHitBoundary(m_hoseWaitUntilHitBoundary);
        p.setWaitUntilFirstHitRBD(m_hoseWaitUntilHitRBD);

        //add new particle to particle list
        m_particleList.push_back(p);
    }

    std::cout << "Hose for fluid : " << m_currentHoseableFluid << "\tnew particle count : " << m_particleList.size() << "\n";

    //reinitialise neighbour structure due to new particles being added
    m_neighbour->reInitialise(m_particleList.size(), m_smoothingLength);
}

void Solver::updateFluid
                (
                    Environment* io_environment,
                    Integration* _integration
                )
{
    //clear neighbour map
    m_neighbour->clearHashmap();

#pragma omp parallel
{
    //refresh neighbour list
    m_neighbour->refreshHashmap(m_particleList);
}

#pragma omp parallel
{
#pragma omp for schedule(dynamic, 50)
    //density-pressure calculations
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        //get neighbours of current particle
        std::vector<FluidParticle> neighbourList = m_neighbour->determineNeighbours(m_particleList[i].getPosition());

        //calculate density and pressure
        ngl::Real density = 0.0;
        for (int j = 0; j < neighbourList.size(); j++)
        {
            //accumulate density
            density += neighbourList[j].getMass() * WKernel_Poly((m_particleList[i].getPosition() - neighbourList[j].getPosition()).length());
        }
        m_particleList[i].setDensity(density);

        //calculate pressure
        m_particleList[i].calculatePressure();
    }

#pragma omp for schedule(dynamic, 50)
    //force calculations
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        //reset net forces accumulate
        m_particleList[i].resetForce();

        //get neighbours of current particle
        std::vector<FluidParticle> neighbourList = m_neighbour->determineNeighbours(m_particleList[i].getPosition());

        //calculate internal forces of current particle
        ngl::Real massPerDensity = 0;
        ngl::Vector separationVector = 0;
        ngl::Real separationLength = 0.0;

        ngl::Vector pressure = 0;
        ngl::Vector viscosity = 0;

        ngl::Vector surfaceGradient = 0;
        ngl::Real surfaceLaplacian = 0;

        ngl::Vector surfaceTensionGradient = 0;
        ngl::Real surfaceTensionLaplacian = 0;
        ngl::Vector interfaceTensionGradient = 0;
        ngl::Real interfaceTensionLaplacian = 0;

        //apply fluid rules if enabled
        if (m_particleList[i].fluidRulesEnabled())
        {
            for (int j = 0; j < neighbourList.size(); j++)
            {
                //get neighbour
                FluidParticle neighbour = neighbourList[j];

                //pre-calculate some repetitive values
                massPerDensity = neighbour.getMass() / neighbour.getDensity();
                separationVector = m_particleList[i].getPosition() - neighbour.getPosition();
                separationLength = separationVector.length();

                //add in pressure and viscosity force only for neighbours, not the current particle itself
                if (neighbour.getId() != m_particleList[i].getId())
                {
                    //accumulate pressure force component
                    pressure = pressure +
                               (
                                       ((m_particleList[i].getPressure() / pow(m_particleList[i].getDensity(), 2)) + (neighbour.getPressure() / pow(neighbour.getDensity(), 2))) *
                                       neighbour.getMass() *
                                       WKernel_Pressure_Gradient(separationVector)
                                       );

                    //accumulate viscous force component
                    viscosity = viscosity +
                                (
                                        /*((currentParticle.getViscosityConstant() + neighbour.getViscosityConstant()) / 2.0) **/
                                        massPerDensity *
                                        (neighbour.getVelocity() - m_particleList[i].getVelocity()) *
                                        WKernel_Viscosity_Laplacian(separationLength)
                                        );
                }

                //accumulate surface tension-interface force component
                surfaceGradient = WKernel_Poly_Gradient(separationVector);
                surfaceLaplacian = WKernel_Poly_Laplacian(separationLength);

                //accumulate surface gradient
                surfaceTensionGradient = surfaceTensionGradient +
                                         (
                                                 massPerDensity *
                                                 m_particleList[i].getSurfaceColorCoefficient() *
                                                 surfaceGradient
                                                 );

                //accumulate surface laplacian
                surfaceTensionLaplacian = surfaceTensionLaplacian +
                                          (
                                                  massPerDensity *
                                                  m_particleList[i].getSurfaceColorCoefficient() *
                                                  surfaceLaplacian
                                                  );

                //accumulate interface gradient
                interfaceTensionGradient = surfaceTensionGradient +
                                           (
                                                   massPerDensity *
                                                   m_particleList[i].getInterfaceColorCoefficient() *
                                                   surfaceGradient
                                                   );

                //accumulate interface laplacian
                interfaceTensionLaplacian = surfaceTensionLaplacian +
                                            (
                                                    massPerDensity *
                                                    m_particleList[i].getInterfaceColorCoefficient() *
                                                    surfaceLaplacian
                                                    );

            }
        }

        //calculate pressure force
        m_particleList[i].setPressureForce(-1.0 * m_particleList[i].getDensity() * pressure);

        //calculate viscosity force
        m_particleList[i].setViscosityForce(viscosity * m_particleList[i].getViscosityConstant());

        //calculate surface tension force
        m_particleList[i].setSurfaceTensionForce(0.0);
        if (surfaceTensionGradient.length() > m_particleList[i].getSurfaceTensionThreshold())
        {
            m_particleList[i].setSurfaceTensionForce(-1.0 * m_particleList[i].getSurfaceTensionCoefficient() * normalise(surfaceTensionGradient) * surfaceTensionLaplacian);
        }

        //calculate interface tension force
        m_particleList[i].setInterfaceTensionForce(0.0);
        if (interfaceTensionGradient.length() > m_particleList[i].getInterfaceTensionThreshold())
        {
            m_particleList[i].setInterfaceTensionForce(-1.0 * m_particleList[i].getInterfaceTensionCoefficient() * normalise(interfaceTensionGradient) * interfaceTensionLaplacian);
        }

        //calculate gravity force (=density * g)
        m_particleList[i].setGravityForce(new ngl::Vector(0, m_particleList[i].getDensity() * Configuration::s_accelerationOfFreeFall, 0));

        //calculate acceleration (= force / mass)
        m_particleList[i].accumulateNetForce(m_particleList[i].getPressureForce());
        m_particleList[i].accumulateNetForce(m_particleList[i].getViscosityForce());
        m_particleList[i].accumulateNetForce(m_particleList[i].getSurfaceTensionForce());
        m_particleList[i].accumulateNetForce(m_particleList[i].getInterfaceTensionForce());
        m_particleList[i].accumulateNetForce(m_particleList[i].getGravityForce());
        m_particleList[i].updateAcceleration(m_particleList[i].getNetForce() / m_particleList[i].getDensity());
    }

#pragma omp for schedule(dynamic, 50)
    //fluid integration and movement
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (m_particleList[i].getMoveable())
        {
            //determine next (position, velocity) of particle
            _integration->integrateNext(m_particleList[i]);
        }
    }

#pragma omp single
    //obstacles movement
    io_environment->updateObstacles(_integration);

#pragma omp for schedule(dynamic, 50)
    //collision check and resolution for fluid particles
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        //check and resolve collision against obstacles and boundary
        bool hitBoundary = false; bool hitRBD = false;
        io_environment->checkAndResolveCollision(m_particleList[i], hitBoundary, hitRBD, false);

        //update first time hit issues
        m_particleList[i].updateHitForTheFirstTime(hitBoundary, hitRBD);
    }
}

//    //increment iteration counter
//    iteration++;
//    std::cout << iteration << "\n";

//    if ((iteration % 1001) == 0)
//    {
//        time_t now; time(&now);
//        printf("%s\n", ctime(&now));
//    }
}

void Solver::render
                (
                    ngl::TransformStack _txStack,
                    ShaderLibrary* io_sman,
                    const std::string _shader
                )
{
    ngl::VBOPrimitives* prim = ngl::VBOPrimitives::instance();

    _txStack.pushTransform();
    {
        for (int i = 0; i < m_particleList.size(); ++i)
        {
            _txStack.getCurrentTransform().setPosition(m_particleList[i].getPosition().m_x, m_particleList[i].getPosition().m_y, m_particleList[i].getPosition().m_z);

            //pass vertex info to shader
            io_sman->updateModel(_shader, _txStack.getCurrAndGlobal().getMatrix(), false);

            //pass color info to shader
            io_sman->updateColor(_shader, m_particleList[i].getColour(), false);

            prim->draw("particleSphere");
        }
    }
    _txStack.popTransform();
}

std::vector<std::vector<ngl::Vector> > Solver::getPositionList()
{
    std::vector<std::vector<ngl::Vector> > positionList;

    //loop through the prototype list that contains 1 particle prototype for each fluid
    for (int i = 0; i < m_hoseParticlePrototypeList.size(); ++i)
    {
        //get name
        std::string name = m_hoseParticlePrototypeList[i].getName();
        std::vector<ngl::Vector> list;
        list.reserve(m_particleList.size());

        //loop through particle list
        for (int j = 0; j < m_particleList.size(); ++j)
        {
            //search for particles having this name
            if (m_particleList[j].getName() == name)
            {
                //add this particle's position to list
                list.push_back(m_particleList[j].getPosition());
            }
        }

        //add this fluid's position list to the big list
        positionList.push_back(list);
    }

    return positionList;
}

std::vector<std::string> Solver::getNameList()
{
    std::vector<std::string> list;
    list.reserve(m_hoseParticlePrototypeList.size());

    for (int i = 0; i < m_hoseParticlePrototypeList.size(); ++i) list.push_back(m_hoseParticlePrototypeList[i].getName());

    return list;
}

void Solver::updateRestDensityForFluid
                (
                    const int _fluidId,
                    const ngl::Real _restDensity
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setMaterialRestDensity(_restDensity);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setMaterialRestDensity(_restDensity);
}

void Solver::updateGasConstantForFluid
                (
                    const int _fluidId,
                    const ngl::Real _gasConstant
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setGasConstant(_gasConstant);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setGasConstant(_gasConstant);
}

void Solver::updateViscosityConstantForFluid
                (
                    const int _fluidId,
                    const ngl::Real _viscosityConstant
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setViscosityConstant(_viscosityConstant);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setViscosityConstant(_viscosityConstant);
}

void Solver::updateSurfaceCoefficientForFluid
                (
                    const int _fluidId,
                    const ngl::Real _coefficient
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setSurfaceTensionCoefficient(_coefficient);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setSurfaceTensionCoefficient(_coefficient);
}

void Solver::updateSurfaceThresholdForFluid
                (
                    const int _fluidId,
                    const ngl::Real _threshold
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setSurfaceTensionThreshold(_threshold);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setSurfaceTensionThreshold(_threshold);
}

void Solver::updateInterfaceCoefficientForFluid
                (
                    const int _fluidId,
                    const ngl::Real _coefficient
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setInterfaceTensionCoefficient(_coefficient);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setInterfaceTensionCoefficient(_coefficient);
}

void Solver::updateInterfaceThresholdForFluid
                (
                    const int _fluidId,
                    const ngl::Real _threshold
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setInterfaceTensionThreshold(_threshold);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setInterfaceTensionThreshold(_threshold);
}

void Solver::updateInterfaceColorForFluid
                (
                    const int _fluidId,
                    const ngl::Real _color
                )
{
    //loop for particles of the concerned fluid
    for (int i = 0; i < m_particleList.size(); ++i)
    {
        if (getFluidIdFromName(m_particleList[i].getName()) == _fluidId)
        {
            //found particles of the concerned fluid -> update it
            m_particleList[i].setInterfaceColorCoefficient(_color);
        }
    }

    m_hoseParticlePrototypeList[_fluidId].setInterfaceColorCoefficient(_color);
}