Precomputed Kernel

[luchete80]

[5] Jan Bender and Dan Koschier. Divergence-free smoothed particle hydrodynamics. In ACM SIGGRAPH / Eurographics Symposium on Computer Animation, SCA '15, pages 147–155, New York, NY, USA, 2015. ACM.

[6] Jan Bender and Dan Koschier. Divergence-free sph for incompressible and viscous fluids. IEEE Transactions on Visualization and Computer Graphics, 23(3):1193–1206, 2017.

[luchete80]

    template<typename KernelType, unsigned int resolution = 10000u>
    class PrecomputedKernel
    {
    protected:
        static Real m_W[resolution];
        static Real m_gradW[resolution + 1];
        static Real m_radius;
        static Real m_radius2;
        static Real m_invStepSize;
        static Real m_W_zero;
    public:
        static Real getRadius() { return m_radius; }
        static void setRadius(Real val)
        {
            m_radius = val;
            m_radius2 = m_radius * m_radius;
            KernelType::setRadius(val);
            const Real stepSize = m_radius / (Real)(resolution-1);
            m_invStepSize = static_cast<Real>(1.0) / stepSize;
            for (unsigned int i = 0; i < resolution; i++)
            {
                const Real posX = stepSize * (Real)i;       // Store kernel values in the middle of an interval
                m_W[i] = KernelType::W(posX);
                KernelType::setRadius(val);
                if (posX > 1.0e-9)
                    m_gradW[i] = KernelType::gradW(Vector3r(posX, 0.0, 0.0))[0] / posX;
                else
                    m_gradW[i] = 0.0;
            }
            m_gradW[resolution] = 0.0;
            m_W_zero = W(static_cast<Real>(0));
        }

    public:
        static Real W(const Vector3r &r)
        {
            Real res = 0.0;
            const Real r2 = r.squaredNorm();
            if (r2 <= m_radius2)
            {
                const Real rl = sqrt(r2);
                const unsigned int pos = std::min<unsigned int>((unsigned int)(rl * m_invStepSize), resolution-2u);
                res = static_cast<Real>(0.5)*(m_W[pos]+ m_W[pos+1]);
            }
            return res;
        }

        static Real W(const Real r)
        {
            Real res = 0.0;
            if (r <= m_radius)
            {
                const unsigned int pos = std::min<unsigned int>((unsigned int)(r * m_invStepSize), resolution-2u);
                res = static_cast<Real>(0.5)*(m_W[pos] + m_W[pos + 1]);
            }
            return res;
        }

        static Vector3r gradW(const Vector3r &r)
        {
            Vector3r res;
            const Real rl = r.norm();
            if (rl <= m_radius)
            {
                //const Real rl = sqrt(r2);
                const unsigned int pos = std::min<unsigned int>(static_cast<unsigned int>(rl * m_invStepSize), resolution-2u);
                res = static_cast<Real>(0.5)*(m_gradW[pos] + m_gradW[pos + 1]) * r;
            }
            else
                res.setZero();

            return res;
        }

        static Real W_zero()
        {
            return m_W_zero;
        }
    };

[luchete80]

Algorithm 1 outlines the simulation with our novel method in detail. First, we determine the particle neighborhoods Ni using compact hashing [Ihmsen et al. 2011]. Then for each particle we compute the density ρi and the factor αi (see Section 3.2). αi is a common factor of both solvers which is required to correct the density error and the divergence error, respectively. Since this factor solely depends on the current positions, it is precomputed before executing the solvers which reduces the computational effort of both solvers significantly.

[luchete80]