MPMSystemData-Impl.hpp

Go to the documentation of this file.

// This code is based on Jet framework.
// Copyright (c) 2018 Doyub Kim
// CubbyFlow is voxel-based fluid simulation engine for computer games.
// Copyright (c) 2020 CubbyFlow Team
// Core Part: Chris Ohk, Junwoo Hwang, Jihong Sin, Seungwoo Yoo
// AI Part: Dongheon Cho, Minseo Kim
// We are making my contributions/submissions to this project solely in our
// personal capacity and are not conveying any rights to any intellectual
// property of any third parties.
 
#ifndef CUBBYFLOW_MPM_SYSTEM_DATA_IMPL_HPP
#define CUBBYFLOW_MPM_SYSTEM_DATA_IMPL_HPP
 
#include <Core/Utils/IterationUtils.hpp>
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <utility>
 
namespace CubbyFlow
{
template <size_t N>
double CubicBSplineKernel<N>::Weight(double x)
{
    const double ax = std::abs(x);
 
    if (ax < 1.0)
    {
        return 0.5 * ax * ax * ax - ax * ax + 2.0 / 3.0;
    }
 
    if (ax < 2.0)
    {
        const double d = 2.0 - ax;
        return d * d * d / 6.0;
    }
 
    return 0.0;
}
 
template <size_t N>
double CubicBSplineKernel<N>::Gradient(double x)
{
    const double ax = std::abs(x);
 
    if (ax < 1.0)
    {
        return x * (1.5 * ax - 2.0);
    }
 
    if (ax < 2.0)
    {
        const double d = 2.0 - ax;
        return -0.5 * d * d * std::copysign(1.0, x);
    }
 
    return 0.0;
}
 
template <size_t N>
void CubicBSplineKernel<N>::GetStencilCoordinates(
    const Vector<double, N>& position, const Vector<double, N>& gridSpacing,
    const Vector<double, N>& dataOrigin, Vector<double, N>* normalized,
    Vector<ssize_t, N>* firstIndex)
{
    for (size_t axis = 0; axis < N; ++axis)
    {
        if (!std::isfinite(position[axis]) ||
            !std::isfinite(gridSpacing[axis]) ||
            !std::isfinite(dataOrigin[axis]) || gridSpacing[axis] <= 0.0)
        {
            throw std::invalid_argument("Invalid cubic B-spline input.");
        }
 
        (*normalized)[axis] =
            (position[axis] - dataOrigin[axis]) / gridSpacing[axis];
 
        const double lowestIndex = std::nextafter(
            static_cast<double>(std::numeric_limits<ssize_t>::lowest()) + 1.0,
            0.0);
        if (const double highestIndex = std::nextafter(
                static_cast<double>(std::numeric_limits<ssize_t>::max()) - 2.0,
                0.0);
            !std::isfinite((*normalized)[axis]) ||
            (*normalized)[axis] < lowestIndex ||
            (*normalized)[axis] > highestIndex)
        {
            throw std::invalid_argument(
                "Cubic B-spline index is out of range.");
        }
 
        (*firstIndex)[axis] =
            static_cast<ssize_t>(std::floor((*normalized)[axis])) - 1;
    }
}
 
template <size_t N>
CubicBSplineKernel<N>::Entry CubicBSplineKernel<N>::GetStencilEntry(
    const Vector<size_t, N>& offset, const Vector<double, N>& normalized,
    const Vector<ssize_t, N>& firstIndex, const Vector<double, N>& gridSpacing)
{
    Entry entry;
    std::array<double, N> axisWeights;
 
    entry.weight = 1.0;
 
    for (size_t axis = 0; axis < N; ++axis)
    {
        entry.index[axis] =
            firstIndex[axis] + static_cast<ssize_t>(offset[axis]);
        axisWeights[axis] =
            Weight(normalized[axis] - static_cast<double>(entry.index[axis]));
        entry.weight *= axisWeights[axis];
    }
 
    for (size_t axis = 0; axis < N; ++axis)
    {
        entry.gradient[axis] =
            Gradient(normalized[axis] -
                     static_cast<double>(entry.index[axis])) /
            gridSpacing[axis];
 
        for (size_t other = 0; other < N; ++other)
        {
            if (other != axis)
            {
                entry.gradient[axis] *= axisWeights[other];
            }
        }
    }
 
    return entry;
}
 
template <size_t N>
CubicBSplineKernel<N>::Stencil CubicBSplineKernel<N>::GetStencil(
    const Vector<double, N>& position, const Vector<double, N>& gridSpacing,
    const Vector<double, N>& dataOrigin)
{
    Vector<double, N> normalized;
    Vector<ssize_t, N> firstIndex;
    GetStencilCoordinates(position, gridSpacing, dataOrigin, &normalized,
                          &firstIndex);
 
    std::array<Entry, STENCIL_SIZE> result;
    size_t flatIndex = 0;
 
    ForEachIndex(Vector<size_t, N>::MakeConstant(4),
                 [&result, &flatIndex, &normalized, &firstIndex,
                  &gridSpacing](auto... rawIndices) {
                     const Vector<size_t, N> offset{ rawIndices... };
                     result[flatIndex++] = GetStencilEntry(
                         offset, normalized, firstIndex, gridSpacing);
                 });
 
    return result;
}
 
template <size_t N>
MPMTransferSystemData<N>::MPMTransferSystemData(
    const Vector<size_t, N>& resolution, const Vector<double, N>& gridSpacing,
    const Vector<double, N>& gridOrigin, size_t numberOfParticles)
    : Base{ numberOfParticles },
      m_particleMasses(numberOfParticles, Base::Mass()),
      m_initialVolumes(numberOfParticles, 0.0)
{
    ResizeGrid(resolution, gridSpacing, gridOrigin);
}
 
template <size_t N>
void MPMTransferSystemData<N>::Resize(size_t newNumberOfParticles)
{
    Base::Resize(newNumberOfParticles);
 
    m_particleMasses.Resize(newNumberOfParticles, Base::Mass());
    m_initialVolumes.Resize(newNumberOfParticles, 0.0);
}
 
template <size_t N>
void MPMTransferSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
{
    Base::Deserialize(buffer);
    ResetTransferState();
}
 
template <size_t N>
void MPMTransferSystemData<N>::Set(const ParticleSystemData<N>& other)
{
    Base::Set(other);
    ResetTransferState();
}
 
template <size_t N>
void MPMTransferSystemData<N>::ResizeGrid(const Vector<size_t, N>& resolution,
                                          const Vector<double, N>& gridSpacing,
                                          const Vector<double, N>& gridOrigin)
{
    ValidateGridParameters(resolution, gridSpacing, gridOrigin);
 
    m_gridMass.Resize(resolution, gridSpacing, gridOrigin);
    m_gridVelocities.Resize(resolution, gridSpacing, gridOrigin);
    m_gridVelocitiesBeforeUpdate.Resize(resolution, gridSpacing, gridOrigin);
}
 
template <size_t N>
ConstArrayView1<double> MPMTransferSystemData<N>::ParticleMasses() const
{
    return m_particleMasses.View();
}
 
template <size_t N>
ArrayView1<double> MPMTransferSystemData<N>::ParticleMasses()
{
    return m_particleMasses.View();
}
 
template <size_t N>
ConstArrayView1<double> MPMTransferSystemData<N>::InitialVolumes() const
{
    return m_initialVolumes.View();
}
 
template <size_t N>
ArrayView1<double> MPMTransferSystemData<N>::InitialVolumes()
{
    return m_initialVolumes.View();
}
 
template <size_t N>
const VertexCenteredScalarGrid<N>& MPMTransferSystemData<N>::GridMass() const
{
    return m_gridMass;
}
 
template <size_t N>
VertexCenteredScalarGrid<N>& MPMTransferSystemData<N>::GridMass()
{
    return m_gridMass;
}
 
template <size_t N>
const VertexCenteredVectorGrid<N>& MPMTransferSystemData<N>::GridVelocities()
    const
{
    return m_gridVelocities;
}
 
template <size_t N>
VertexCenteredVectorGrid<N>& MPMTransferSystemData<N>::GridVelocities()
{
    return m_gridVelocities;
}
 
template <size_t N>
const VertexCenteredVectorGrid<N>&
MPMTransferSystemData<N>::GridVelocitiesBeforeUpdate() const
{
    return m_gridVelocitiesBeforeUpdate;
}
 
template <size_t N>
VertexCenteredVectorGrid<N>&
MPMTransferSystemData<N>::GridVelocitiesBeforeUpdate()
{
    return m_gridVelocitiesBeforeUpdate;
}
 
template <size_t N>
double MPMTransferSystemData<N>::FLIPBlendingFactor() const
{
    return m_flipBlendingFactor;
}
 
template <size_t N>
void MPMTransferSystemData<N>::SetFLIPBlendingFactor(double factor)
{
    if (!std::isfinite(factor) || factor < 0.0 || factor > 1.0)
    {
        throw std::invalid_argument("FLIP blending factor must be in [0, 1].");
    }
 
    m_flipBlendingFactor = factor;
}
 
template <size_t N>
void MPMTransferSystemData<N>::TransferFromParticlesToGrid()
{
    ValidateGridState();
 
    const auto positions = this->Positions();
    const auto velocities = this->Velocities();
 
    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
    {
        if (!std::isfinite(m_particleMasses[i]) || m_particleMasses[i] <= 0.0 ||
            !IsFinite(positions[i]) || !IsFinite(velocities[i]))
        {
            throw std::invalid_argument("Invalid MPM particle state.");
        }
    }
 
    VertexCenteredScalarGrid<N> nextGridMass{ m_gridMass.Resolution(),
                                              m_gridMass.GridSpacing(),
                                              m_gridMass.Origin() };
    VertexCenteredVectorGrid<N> nextGridVelocities{
        m_gridVelocities.Resolution(), m_gridVelocities.GridSpacing(),
        m_gridVelocities.Origin()
    };
 
    const auto dataSize = nextGridMass.DataSize();
    const auto gridSpacing = nextGridMass.GridSpacing();
    const auto dataOrigin = nextGridMass.DataOrigin();
 
    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
    {
        const auto stencil = CubicBSplineKernel<N>::GetStencil(
            positions[i], gridSpacing, dataOrigin);
 
        for (const auto& entry : stencil)
        {
            if (entry.weight == 0.0)
            {
                continue;
            }
 
            const auto index = ClampIndex(entry.index, dataSize);
            const double mass = entry.weight * m_particleMasses[i];
 
            nextGridMass(index) += mass;
            nextGridVelocities(index) += mass * velocities[i];
        }
    }
 
    nextGridMass.ForEachDataPointIndex([&nextGridMass, &nextGridVelocities](
                                           const Vector<size_t, N>& index) {
        const double mass = nextGridMass(index);
 
        if (!std::isfinite(mass) ||
            !MPMTransferSystemData<N>::IsFinite(nextGridVelocities(index)))
        {
            throw std::invalid_argument("Invalid MPM grid accumulation.");
        }
 
        if (mass > 0.0)
        {
            nextGridVelocities(index) /= mass;
 
            if (!MPMTransferSystemData<N>::IsFinite(nextGridVelocities(index)))
            {
                throw std::invalid_argument(
                    "Invalid MPM grid velocity update.");
            }
        }
    });
 
    VertexCenteredVectorGrid<N> nextGridVelocitiesBeforeUpdate =
        nextGridVelocities;
 
    m_gridMass = std::move(nextGridMass);
    m_gridVelocities = std::move(nextGridVelocities);
    m_gridVelocitiesBeforeUpdate = std::move(nextGridVelocitiesBeforeUpdate);
}
 
template <size_t N>
void MPMTransferSystemData<N>::TransferFromGridToParticles()
{
    ValidateGridToParticleState();
    TransferFromGridToParticlesUnchecked();
}
 
template <size_t N>
void MPMTransferSystemData<N>::ValidateGridToParticleState() const
{
    ValidateGridState();
 
    const auto positions = this->Positions();
    const auto velocities = this->Velocities();
 
    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
    {
        if (!IsFinite(positions[i]) || !IsFinite(velocities[i]))
        {
            throw std::invalid_argument("Invalid MPM particle state.");
        }
    }
 
    m_gridVelocities.ForEachDataPointIndex(
        [this](const Vector<size_t, N>& index) {
            if (!IsFinite(m_gridVelocities(index)) ||
                !IsFinite(m_gridVelocitiesBeforeUpdate(index)))
            {
                throw std::invalid_argument("Invalid MPM grid velocity.");
            }
        });
}
 
template <size_t N>
void MPMTransferSystemData<N>::TransferFromGridToParticlesUnchecked()
{
    const auto positions = this->Positions();
    auto velocities = this->Velocities();
    Array1<Vector<double, N>> nextVelocities(this->NumberOfParticles());
 
    const auto dataSize = m_gridVelocities.DataSize();
    const auto gridSpacing = m_gridVelocities.GridSpacing();
    const auto dataOrigin = m_gridVelocities.DataOrigin();
 
    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
    {
        Vector<double, N> picVelocity;
        Vector<double, N> flipDelta;
        const auto stencil = CubicBSplineKernel<N>::GetStencil(
            positions[i], gridSpacing, dataOrigin);
 
        for (const auto& entry : stencil)
        {
            if (entry.weight == 0.0)
            {
                continue;
            }
 
            const auto index = ClampIndex(entry.index, dataSize);
            picVelocity += entry.weight * m_gridVelocities(index);
            flipDelta += entry.weight * (m_gridVelocities(index) -
                                         m_gridVelocitiesBeforeUpdate(index));
        }
 
        const Vector<double, N> flipVelocity = velocities[i] + flipDelta;
        const Vector<double, N> result =
            (1.0 - m_flipBlendingFactor) * picVelocity +
            m_flipBlendingFactor * flipVelocity;
 
        if (!IsFinite(result))
        {
            throw std::invalid_argument(
                "Invalid MPM particle velocity update.");
        }
 
        nextVelocities[i] = result;
    }
 
    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
    {
        velocities[i] = nextVelocities[i];
    }
}
 
template <size_t N>
Matrix<double, N, N> MPMTransferSystemData<N>::ComputeVelocityGradient(
    size_t particleIndex) const
{
    const auto positions = this->Positions();
    const auto dataSize = m_gridVelocities.DataSize();
    const auto stencil = CubicBSplineKernel<N>::GetStencil(
        positions[particleIndex], m_gridVelocities.GridSpacing(),
        m_gridVelocities.DataOrigin());
    Matrix<double, N, N> result;
 
    for (const auto& entry : stencil)
    {
        if (entry.weight == 0.0)
        {
            continue;
        }
 
        const Vector<double, N> velocity =
            m_gridVelocities(ClampIndex(entry.index, dataSize));
 
        for (size_t row = 0; row < N; ++row)
        {
            for (size_t column = 0; column < N; ++column)
            {
                result(row, column) += velocity[row] * entry.gradient[column];
            }
        }
    }
 
    return result;
}
 
template <size_t N>
Vector<size_t, N> MPMTransferSystemData<N>::ClampIndex(
    const Vector<ssize_t, N>& index, const Vector<size_t, N>& dataSize)
{
    Vector<size_t, N> result;
 
    for (size_t axis = 0; axis < N; ++axis)
    {
        result[axis] = index[axis] < 0
                           ? 0
                           : std::min(static_cast<size_t>(index[axis]),
                                      dataSize[axis] - 1);
    }
 
    return result;
}
 
template <size_t N>
bool MPMTransferSystemData<N>::IsFinite(const Vector<double, N>& value)
{
    for (size_t axis = 0; axis < N; ++axis)
    {
        if (!std::isfinite(value[axis]))
        {
            return false;
        }
    }
 
    return true;
}
 
template <size_t N>
void MPMTransferSystemData<N>::ValidateGridParameters(
    const Vector<size_t, N>& resolution, const Vector<double, N>& gridSpacing,
    const Vector<double, N>& gridOrigin)
{
    for (size_t axis = 0; axis < N; ++axis)
    {
        if (resolution[axis] == 0 ||
            resolution[axis] == std::numeric_limits<size_t>::max() ||
            !std::isfinite(gridSpacing[axis]) || gridSpacing[axis] <= 0.0 ||
            !std::isfinite(gridOrigin[axis]))
        {
            throw std::invalid_argument("Invalid MPM grid parameters.");
        }
    }
}
 
template <size_t N>
void MPMTransferSystemData<N>::ValidateGridState() const
{
    ValidateGridParameters(m_gridMass.Resolution(), m_gridMass.GridSpacing(),
                           m_gridMass.Origin());
    ValidateGridParameters(m_gridVelocities.Resolution(),
                           m_gridVelocities.GridSpacing(),
                           m_gridVelocities.Origin());
    ValidateGridParameters(m_gridVelocitiesBeforeUpdate.Resolution(),
                           m_gridVelocitiesBeforeUpdate.GridSpacing(),
                           m_gridVelocitiesBeforeUpdate.Origin());
 
    if (m_gridMass.Resolution() != m_gridVelocities.Resolution() ||
        m_gridMass.Resolution() != m_gridVelocitiesBeforeUpdate.Resolution() ||
        m_gridMass.GridSpacing() != m_gridVelocities.GridSpacing() ||
        m_gridMass.GridSpacing() !=
            m_gridVelocitiesBeforeUpdate.GridSpacing() ||
        m_gridMass.Origin() != m_gridVelocities.Origin() ||
        m_gridMass.Origin() != m_gridVelocitiesBeforeUpdate.Origin())
    {
        throw std::invalid_argument("MPM grids must have matching geometry.");
    }
}
 
template <size_t N>
void MPMTransferSystemData<N>::ResetTransferState()
{
    m_particleMasses.Fill(Base::Mass());
    m_initialVolumes.Fill(0.0);
    m_gridMass.Fill(0.0, ExecutionPolicy::Serial);
    m_gridVelocities.Fill(Vector<double, N>{}, ExecutionPolicy::Serial);
    m_gridVelocitiesBeforeUpdate.Fill(Vector<double, N>{},
                                      ExecutionPolicy::Serial);
}
 
template <size_t N>
MPMSystemData<N>::MPMSystemData(const Vector<size_t, N>& resolution,
                                const Vector<double, N>& gridSpacing,
                                const Vector<double, N>& gridOrigin,
                                size_t numberOfParticles)
    : TransferBase{ resolution, gridSpacing, gridOrigin, numberOfParticles },
      m_deformationStates{ numberOfParticles, DeformationState{} }
{
    // Do nothing
}
 
template <size_t N>
void MPMSystemData<N>::Resize(size_t newNumberOfParticles)
{
    TransferBase::Resize(newNumberOfParticles);
    m_deformationStates.Resize(newNumberOfParticles, DeformationState{});
}
 
template <size_t N>
void MPMSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
{
    TransferBase::Deserialize(buffer);
    m_deformationStates.Fill(DeformationState{});
}
 
template <size_t N>
void MPMSystemData<N>::Set(const ParticleSystemData<N>& other)
{
    TransferBase::Set(other);
    m_deformationStates.Fill(DeformationState{});
}
 
template <size_t N>
ConstArrayView1<typename MPMSystemData<N>::DeformationState>
MPMSystemData<N>::DeformationStates() const
{
    return m_deformationStates.View();
}
 
template <size_t N>
ArrayView1<typename MPMSystemData<N>::DeformationState>
MPMSystemData<N>::DeformationStates()
{
    return m_deformationStates.View();
}
}  // namespace CubbyFlow
 
#endif