From d0040a8ce12c62f16415164d6e5859cfbe5475a2 Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Mon, 24 Aug 2026 21:05:07 +0900
Subject: [PATCH 1/8] refactor: extract shared MPM transfer state

---
 .../Core/Particle/MPM/MPMSystemData-Impl.hpp  | 241 +++++++++++++-----
 Includes/Core/Particle/MPM/MPMSystemData.hpp  |  84 ++++--
 2 files changed, 247 insertions(+), 78 deletions(-)

diff --git a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
index 0577661083..3232811c66 100644
--- a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
+++ b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
@@ -17,6 +17,7 @@
 #include <cmath>
 #include <limits>
 #include <stdexcept>
+#include <utility>
 
 namespace CubbyFlow
 {
@@ -158,10 +159,9 @@ CubicBSplineKernel<N>::Stencil CubicBSplineKernel<N>::GetStencil(
 }
 
 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)
+MPMTransferSystemData<N>::MPMTransferSystemData(
+    const Vector<size_t, N>& resolution, const Vector<double, N>& gridSpacing,
+    const Vector<double, N>& gridOrigin, size_t numberOfParticles)
     : Base{}
 {
     ResizeGrid(resolution, gridSpacing, gridOrigin);
@@ -169,33 +169,32 @@ MPMSystemData<N>::MPMSystemData(const Vector<size_t, N>& resolution,
 }
 
 template <size_t N>
-void MPMSystemData<N>::Resize(size_t newNumberOfParticles)
+void MPMTransferSystemData<N>::Resize(size_t newNumberOfParticles)
 {
     Base::Resize(newNumberOfParticles);
 
     m_particleMasses.Resize(newNumberOfParticles, Base::Mass());
     m_initialVolumes.Resize(newNumberOfParticles, 0.0);
-    m_deformationStates.Resize(newNumberOfParticles, DeformationState{});
 }
 
 template <size_t N>
-void MPMSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
+void MPMTransferSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
 {
     Base::Deserialize(buffer);
-    ResetMPMState();
+    ResetTransferState();
 }
 
 template <size_t N>
-void MPMSystemData<N>::Set(const ParticleSystemData<N>& other)
+void MPMTransferSystemData<N>::Set(const ParticleSystemData<N>& other)
 {
     Base::Set(other);
-    ResetMPMState();
+    ResetTransferState();
 }
 
 template <size_t N>
-void MPMSystemData<N>::ResizeGrid(const Vector<size_t, N>& resolution,
-                                  const Vector<double, N>& gridSpacing,
-                                  const Vector<double, N>& gridOrigin)
+void MPMTransferSystemData<N>::ResizeGrid(const Vector<size_t, N>& resolution,
+                                          const Vector<double, N>& gridSpacing,
+                                          const Vector<double, N>& gridOrigin)
 {
     ValidateGridParameters(resolution, gridSpacing, gridOrigin);
 
@@ -205,88 +204,76 @@ void MPMSystemData<N>::ResizeGrid(const Vector<size_t, N>& resolution,
 }
 
 template <size_t N>
-ConstArrayView1<double> MPMSystemData<N>::ParticleMasses() const
+ConstArrayView1<double> MPMTransferSystemData<N>::ParticleMasses() const
 {
     return m_particleMasses.View();
 }
 
 template <size_t N>
-ArrayView1<double> MPMSystemData<N>::ParticleMasses()
+ArrayView1<double> MPMTransferSystemData<N>::ParticleMasses()
 {
     return m_particleMasses.View();
 }
 
 template <size_t N>
-ConstArrayView1<double> MPMSystemData<N>::InitialVolumes() const
+ConstArrayView1<double> MPMTransferSystemData<N>::InitialVolumes() const
 {
     return m_initialVolumes.View();
 }
 
 template <size_t N>
-ArrayView1<double> MPMSystemData<N>::InitialVolumes()
+ArrayView1<double> MPMTransferSystemData<N>::InitialVolumes()
 {
     return m_initialVolumes.View();
 }
 
 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();
-}
-
-template <size_t N>
-const VertexCenteredScalarGrid<N>& MPMSystemData<N>::GridMass() const
+const VertexCenteredScalarGrid<N>& MPMTransferSystemData<N>::GridMass() const
 {
     return m_gridMass;
 }
 
 template <size_t N>
-VertexCenteredScalarGrid<N>& MPMSystemData<N>::GridMass()
+VertexCenteredScalarGrid<N>& MPMTransferSystemData<N>::GridMass()
 {
     return m_gridMass;
 }
 
 template <size_t N>
-const VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocities() const
+const VertexCenteredVectorGrid<N>& MPMTransferSystemData<N>::GridVelocities()
+    const
 {
     return m_gridVelocities;
 }
 
 template <size_t N>
-VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocities()
+VertexCenteredVectorGrid<N>& MPMTransferSystemData<N>::GridVelocities()
 {
     return m_gridVelocities;
 }
 
 template <size_t N>
 const VertexCenteredVectorGrid<N>&
-MPMSystemData<N>::GridVelocitiesBeforeUpdate() const
+MPMTransferSystemData<N>::GridVelocitiesBeforeUpdate() const
 {
     return m_gridVelocitiesBeforeUpdate;
 }
 
 template <size_t N>
-VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocitiesBeforeUpdate()
+VertexCenteredVectorGrid<N>&
+MPMTransferSystemData<N>::GridVelocitiesBeforeUpdate()
 {
     return m_gridVelocitiesBeforeUpdate;
 }
 
 template <size_t N>
-double MPMSystemData<N>::FLIPBlendingFactor() const
+double MPMTransferSystemData<N>::FLIPBlendingFactor() const
 {
     return m_flipBlendingFactor;
 }
 
 template <size_t N>
-void MPMSystemData<N>::SetFLIPBlendingFactor(double factor)
+void MPMTransferSystemData<N>::SetFLIPBlendingFactor(double factor)
 {
     if (!std::isfinite(factor) || factor < 0.0 || factor > 1.0)
     {
@@ -297,7 +284,7 @@ void MPMSystemData<N>::SetFLIPBlendingFactor(double factor)
 }
 
 template <size_t N>
-void MPMSystemData<N>::TransferFromParticlesToGrid()
+void MPMTransferSystemData<N>::TransferFromParticlesToGrid()
 {
     ValidateGridState();
 
@@ -313,12 +300,17 @@ void MPMSystemData<N>::TransferFromParticlesToGrid()
         }
     }
 
-    m_gridMass.Fill(0.0, ExecutionPolicy::Serial);
-    m_gridVelocities.Fill(Vector<double, N>{}, ExecutionPolicy::Serial);
+    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 = m_gridMass.DataSize();
-    const auto gridSpacing = m_gridMass.GridSpacing();
-    const auto dataOrigin = m_gridMass.DataOrigin();
+    const auto dataSize = nextGridMass.DataSize();
+    const auto gridSpacing = nextGridMass.GridSpacing();
+    const auto dataOrigin = nextGridMass.DataOrigin();
 
     for (size_t i = 0; i < this->NumberOfParticles(); ++i)
     {
@@ -335,29 +327,55 @@ void MPMSystemData<N>::TransferFromParticlesToGrid()
             const auto index = ClampIndex(entry.index, dataSize);
             const double mass = entry.weight * m_particleMasses[i];
 
-            m_gridMass(index) += mass;
-            m_gridVelocities(index) += mass * velocities[i];
+            nextGridMass(index) += mass;
+            nextGridVelocities(index) += mass * velocities[i];
         }
     }
 
-    m_gridMass.ForEachDataPointIndex([this](const Vector<size_t, N>& index) {
-        const double mass = m_gridMass(index);
+    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)
         {
-            m_gridVelocities(index) /= mass;
+            nextGridVelocities(index) /= mass;
+
+            if (!MPMTransferSystemData<N>::IsFinite(nextGridVelocities(index)))
+            {
+                throw std::invalid_argument(
+                    "Invalid MPM grid velocity update.");
+            }
         }
     });
 
-    m_gridVelocitiesBeforeUpdate.Set(m_gridVelocities);
+    VertexCenteredVectorGrid<N> nextGridVelocitiesBeforeUpdate =
+        nextGridVelocities;
+
+    m_gridMass = std::move(nextGridMass);
+    m_gridVelocities = std::move(nextGridVelocities);
+    m_gridVelocitiesBeforeUpdate = std::move(nextGridVelocitiesBeforeUpdate);
 }
 
 template <size_t N>
-void MPMSystemData<N>::TransferFromGridToParticles()
+void MPMTransferSystemData<N>::TransferFromGridToParticles()
+{
+    ValidateGridToParticleState();
+    TransferFromGridToParticlesUnchecked();
+}
+
+template <size_t N>
+void MPMTransferSystemData<N>::ValidateGridToParticleState() const
 {
     ValidateGridState();
 
     const auto positions = this->Positions();
-    auto velocities = this->Velocities();
+    const auto velocities = this->Velocities();
 
     for (size_t i = 0; i < this->NumberOfParticles(); ++i)
     {
@@ -375,6 +393,14 @@ void MPMSystemData<N>::TransferFromGridToParticles()
                 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();
@@ -405,12 +431,56 @@ void MPMSystemData<N>::TransferFromGridToParticles()
             (1.0 - m_flipBlendingFactor) * picVelocity +
             m_flipBlendingFactor * flipVelocity;
 
-        velocities[i] = result;
+        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>
-Vector<size_t, N> MPMSystemData<N>::ClampIndex(
+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;
@@ -427,7 +497,7 @@ Vector<size_t, N> MPMSystemData<N>::ClampIndex(
 }
 
 template <size_t N>
-bool MPMSystemData<N>::IsFinite(const Vector<double, N>& value)
+bool MPMTransferSystemData<N>::IsFinite(const Vector<double, N>& value)
 {
     for (size_t axis = 0; axis < N; ++axis)
     {
@@ -441,7 +511,7 @@ bool MPMSystemData<N>::IsFinite(const Vector<double, N>& value)
 }
 
 template <size_t N>
-void MPMSystemData<N>::ValidateGridParameters(
+void MPMTransferSystemData<N>::ValidateGridParameters(
     const Vector<size_t, N>& resolution, const Vector<double, N>& gridSpacing,
     const Vector<double, N>& gridOrigin)
 {
@@ -458,7 +528,7 @@ void MPMSystemData<N>::ValidateGridParameters(
 }
 
 template <size_t N>
-void MPMSystemData<N>::ValidateGridState() const
+void MPMTransferSystemData<N>::ValidateGridState() const
 {
     ValidateGridParameters(m_gridMass.Resolution(), m_gridMass.GridSpacing(),
                            m_gridMass.Origin());
@@ -482,16 +552,67 @@ void MPMSystemData<N>::ValidateGridState() const
 }
 
 template <size_t N>
-void MPMSystemData<N>::ResetMPMState()
+void MPMTransferSystemData<N>::ResetTransferState()
 {
     m_particleMasses.Fill(Base::Mass());
     m_initialVolumes.Fill(0.0);
-    m_deformationStates.Fill(DeformationState{});
     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>
+void MPMSystemData<N>::TransferFromGridToParticles()
+{
+    TransferBase::TransferFromGridToParticles();
+}
+
+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
diff --git a/Includes/Core/Particle/MPM/MPMSystemData.hpp b/Includes/Core/Particle/MPM/MPMSystemData.hpp
index 008c6ca6bf..25b5c4407c 100644
--- a/Includes/Core/Particle/MPM/MPMSystemData.hpp
+++ b/Includes/Core/Particle/MPM/MPMSystemData.hpp
@@ -62,22 +62,21 @@ class CubicBSplineKernel final
 };
 
 //!
-//! \brief N-D material point method particle and grid state.
+//! \brief Shared N-D material point method particle-grid transfer state.
 //!
 //! \note Serialization preserves only the inherited particle-system state;
-//! MPM-specific particle and grid state is reset after deserialization.
+//! transfer-specific particle and grid state is reset after deserialization.
 //!
 template <size_t N>
-class MPMSystemData final : public ParticleSystemData<N>
+class MPMTransferSystemData : public ParticleSystemData<N>
 {
  public:
     static_assert(N == 2 || N == 3, "MPM supports only 2-D and 3-D.");
 
     using Base = ParticleSystemData<N>;
-    using DeformationState = SnowDeformationState<N>;
 
-    //! Constructs MPM state with a vertex-centered background grid.
-    explicit MPMSystemData(
+    //! Constructs transfer state with a vertex-centered background grid.
+    explicit MPMTransferSystemData(
         const Vector<size_t, N>& resolution =
             Vector<size_t, N>::MakeConstant(1),
         const Vector<double, N>& gridSpacing =
@@ -85,13 +84,13 @@ class MPMSystemData final : public ParticleSystemData<N>
         const Vector<double, N>& gridOrigin = Vector<double, N>{},
         size_t numberOfParticles = 0);
 
-    //! Resizes particle state, initializing new MPM attributes.
+    //! Resizes particle state, initializing new transfer attributes.
     void Resize(size_t newNumberOfParticles) override;
 
-    //! Deserializes inherited particle state and resets MPM-specific state.
+    //! Deserializes inherited particle state and resets transfer state.
     void Deserialize(const std::vector<uint8_t>& buffer) override;
 
-    //! Copies inherited particle state and resets MPM-specific state.
+    //! Copies inherited particle state and resets transfer state.
     void Set(const ParticleSystemData<N>& other) override;
 
     //! Resizes the background grid without changing particle state.
@@ -111,12 +110,6 @@ class MPMSystemData final : public ParticleSystemData<N>
     //! Returns per-particle initial volumes.
     [[nodiscard]] ArrayView1<double> InitialVolumes();
 
-    //! Returns per-particle deformation states.
-    [[nodiscard]] ConstArrayView1<DeformationState> DeformationStates() const;
-
-    //! Returns per-particle deformation states.
-    [[nodiscard]] ArrayView1<DeformationState> DeformationStates();
-
     //! Returns grid mass.
     [[nodiscard]] const VertexCenteredScalarGrid<N>& GridMass() const;
 
@@ -146,11 +139,11 @@ class MPMSystemData final : public ParticleSystemData<N>
     //! Stencil nodes outside the finite grid are clamped to its boundary.
     void TransferFromParticlesToGrid();
 
+ protected:
     //! Transfers updated background-grid velocities to particles.
     //! Stencil nodes outside the finite grid are clamped to its boundary.
     void TransferFromGridToParticles();
 
- private:
     [[nodiscard]] static Vector<size_t, N> ClampIndex(
         const Vector<ssize_t, N>& index, const Vector<size_t, N>& dataSize);
 
@@ -162,17 +155,72 @@ class MPMSystemData final : public ParticleSystemData<N>
 
     void ValidateGridState() const;
 
-    void ResetMPMState();
+    void ValidateGridToParticleState() const;
+
+    void TransferFromGridToParticlesUnchecked();
+
+    [[nodiscard]] Matrix<double, N, N> ComputeVelocityGradient(
+        size_t particleIndex) const;
+
+ private:
+    void ResetTransferState();
 
     Array1<double> m_particleMasses;
     Array1<double> m_initialVolumes;
-    Array1<DeformationState> m_deformationStates;
     VertexCenteredScalarGrid<N> m_gridMass;
     VertexCenteredVectorGrid<N> m_gridVelocities;
     VertexCenteredVectorGrid<N> m_gridVelocitiesBeforeUpdate;
     double m_flipBlendingFactor = 0.95;
 };
 
+//!
+//! \brief N-D snow material point method particle and grid state.
+//!
+//! \note Serialization preserves only the inherited particle-system state;
+//! MPM-specific particle and grid state is reset after deserialization.
+//!
+template <size_t N>
+class MPMSystemData final : public MPMTransferSystemData<N>
+{
+ public:
+    static_assert(N == 2 || N == 3, "MPM supports only 2-D and 3-D.");
+
+    using Base = ParticleSystemData<N>;
+    using TransferBase = MPMTransferSystemData<N>;
+    using DeformationState = SnowDeformationState<N>;
+
+    //! Constructs snow MPM state with a vertex-centered background grid.
+    explicit MPMSystemData(
+        const Vector<size_t, N>& resolution =
+            Vector<size_t, N>::MakeConstant(1),
+        const Vector<double, N>& gridSpacing =
+            Vector<double, N>::MakeConstant(1.0),
+        const Vector<double, N>& gridOrigin = Vector<double, N>{},
+        size_t numberOfParticles = 0);
+
+    //! Resizes particle state, initializing new deformation states.
+    void Resize(size_t newNumberOfParticles) override;
+
+    //! Deserializes inherited particle state and resets deformation state.
+    void Deserialize(const std::vector<uint8_t>& buffer) override;
+
+    //! Copies inherited particle state and resets deformation state.
+    void Set(const ParticleSystemData<N>& other) override;
+
+    //! Transfers updated background-grid velocities to particles.
+    //! Stencil nodes outside the finite grid are clamped to its boundary.
+    void TransferFromGridToParticles();
+
+    //! Returns per-particle deformation states.
+    [[nodiscard]] ConstArrayView1<DeformationState> DeformationStates() const;
+
+    //! Returns per-particle deformation states.
+    [[nodiscard]] ArrayView1<DeformationState> DeformationStates();
+
+ private:
+    Array1<DeformationState> m_deformationStates;
+};
+
 //! 2-D material point method system data.
 using MPMSystemData2 = MPMSystemData<2>;
 

From 5693a11eb13def4738f6532d7ae308c54737a5ea Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Mon, 24 Aug 2026 21:05:48 +0900
Subject: [PATCH 2/8] feat: add fluid MPM transfer state

---
 .../Particle/MPM/MPMFluidSystemData-Impl.hpp  | 133 +++++
 .../Core/Particle/MPM/MPMFluidSystemData.hpp  |  93 ++++
 Tests/UnitTests/MPMFluidSystemDataTests.cpp   | 474 ++++++++++++++++++
 3 files changed, 700 insertions(+)
 create mode 100644 Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
 create mode 100644 Includes/Core/Particle/MPM/MPMFluidSystemData.hpp
 create mode 100644 Tests/UnitTests/MPMFluidSystemDataTests.cpp

diff --git a/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
new file mode 100644
index 0000000000..f8787b746e
--- /dev/null
+++ b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
@@ -0,0 +1,133 @@
+// 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_FLUID_SYSTEM_DATA_IMPL_HPP
+#define CUBBYFLOW_MPM_FLUID_SYSTEM_DATA_IMPL_HPP
+
+#include <cmath>
+#include <stdexcept>
+
+namespace CubbyFlow
+{
+template <size_t N>
+MPMFluidSystemData<N>::MPMFluidSystemData(const Vector<size_t, N>& resolution,
+                                          const Vector<double, N>& gridSpacing,
+                                          const Vector<double, N>& gridOrigin,
+                                          size_t numberOfParticles)
+    : Base{ resolution, gridSpacing, gridOrigin, numberOfParticles },
+      m_volumeRatios(numberOfParticles, 1.0),
+      m_velocityGradients(numberOfParticles, MatrixType{})
+{
+    // Do nothing
+}
+
+template <size_t N>
+void MPMFluidSystemData<N>::Resize(size_t newNumberOfParticles)
+{
+    Base::Resize(newNumberOfParticles);
+
+    m_volumeRatios.Resize(newNumberOfParticles, 1.0);
+    m_velocityGradients.Resize(newNumberOfParticles, MatrixType{});
+}
+
+template <size_t N>
+void MPMFluidSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
+{
+    Base::Deserialize(buffer);
+    ResetFluidState();
+}
+
+template <size_t N>
+void MPMFluidSystemData<N>::Set(const ParticleSystemData<N>& other)
+{
+    Base::Set(other);
+    ResetFluidState();
+}
+
+template <size_t N>
+ConstArrayView1<double> MPMFluidSystemData<N>::VolumeRatios() const
+{
+    return m_volumeRatios.View();
+}
+
+template <size_t N>
+ArrayView1<double> MPMFluidSystemData<N>::VolumeRatios()
+{
+    return m_volumeRatios.View();
+}
+
+template <size_t N>
+ConstArrayView1<typename MPMFluidSystemData<N>::MatrixType>
+MPMFluidSystemData<N>::VelocityGradients() const
+{
+    return m_velocityGradients.View();
+}
+
+template <size_t N>
+ArrayView1<typename MPMFluidSystemData<N>::MatrixType>
+MPMFluidSystemData<N>::VelocityGradients()
+{
+    return m_velocityGradients.View();
+}
+
+template <size_t N>
+void MPMFluidSystemData<N>::TransferFromGridToParticles(
+    double timeStepInSeconds)
+{
+    if (!std::isfinite(timeStepInSeconds) || timeStepInSeconds < 0.0)
+    {
+        throw std::invalid_argument("Invalid MPM fluid time step.");
+    }
+
+    this->ValidateGridToParticleState();
+
+    Array1<double> nextVolumeRatios(this->NumberOfParticles(), 1.0);
+    Array1<MatrixType> nextVelocityGradients(this->NumberOfParticles(),
+                                             MatrixType{});
+
+    for (size_t i = 0; i < this->NumberOfParticles(); ++i)
+    {
+        if (!std::isfinite(m_volumeRatios[i]) || m_volumeRatios[i] <= 0.0)
+        {
+            throw std::invalid_argument("Invalid MPM fluid volume ratio.");
+        }
+
+        const MatrixType gradient = this->ComputeVelocityGradient(i);
+        const double volumeIncrement =
+            (MatrixType::MakeIdentity() + timeStepInSeconds * gradient)
+                .Determinant();
+        const double volumeRatio = m_volumeRatios[i] * volumeIncrement;
+
+        if (!std::isfinite(volumeIncrement) || volumeIncrement <= 0.0 ||
+            !std::isfinite(volumeRatio) || volumeRatio <= 0.0)
+        {
+            throw std::invalid_argument("Invalid MPM fluid volume update.");
+        }
+
+        nextVelocityGradients[i] = gradient;
+        nextVolumeRatios[i] = volumeRatio;
+    }
+
+    this->TransferFromGridToParticlesUnchecked();
+
+    m_volumeRatios.Swap(nextVolumeRatios);
+    m_velocityGradients.Swap(nextVelocityGradients);
+}
+
+template <size_t N>
+void MPMFluidSystemData<N>::ResetFluidState()
+{
+    m_volumeRatios.Fill(1.0);
+    m_velocityGradients.Fill(MatrixType{});
+}
+
+}  // namespace CubbyFlow
+
+#endif
diff --git a/Includes/Core/Particle/MPM/MPMFluidSystemData.hpp b/Includes/Core/Particle/MPM/MPMFluidSystemData.hpp
new file mode 100644
index 0000000000..fac13d3d44
--- /dev/null
+++ b/Includes/Core/Particle/MPM/MPMFluidSystemData.hpp
@@ -0,0 +1,93 @@
+// 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_FLUID_SYSTEM_DATA_HPP
+#define CUBBYFLOW_MPM_FLUID_SYSTEM_DATA_HPP
+
+#include <Core/Particle/MPM/MPMSystemData.hpp>
+
+namespace CubbyFlow
+{
+//!
+//! \brief N-D weakly compressible MPM fluid transfer state.
+//!
+//! Stores the particle volume ratio and velocity gradient required by a fluid
+//! constitutive model without carrying snow deformation or plasticity state.
+//!
+template <size_t N>
+class MPMFluidSystemData final : public MPMTransferSystemData<N>
+{
+ public:
+    static_assert(N == 2 || N == 3, "MPM supports only 2-D and 3-D.");
+
+    using Base = MPMTransferSystemData<N>;
+    using MatrixType = Matrix<double, N, N>;
+
+    //! Constructs fluid MPM state with a vertex-centered background grid.
+    explicit MPMFluidSystemData(
+        const Vector<size_t, N>& resolution =
+            Vector<size_t, N>::MakeConstant(1),
+        const Vector<double, N>& gridSpacing =
+            Vector<double, N>::MakeConstant(1.0),
+        const Vector<double, N>& gridOrigin = Vector<double, N>{},
+        size_t numberOfParticles = 0);
+
+    //! Resizes particle state, initializing new fluid attributes.
+    void Resize(size_t newNumberOfParticles) override;
+
+    //! Deserializes inherited particle state and resets fluid state.
+    void Deserialize(const std::vector<uint8_t>& buffer) override;
+
+    //! Copies inherited particle state and resets fluid state.
+    void Set(const ParticleSystemData<N>& other) override;
+
+    //! Returns current-to-reference particle volume ratios.
+    [[nodiscard]] ConstArrayView1<double> VolumeRatios() const;
+
+    //! Returns current-to-reference particle volume ratios.
+    [[nodiscard]] ArrayView1<double> VolumeRatios();
+
+    //! Returns particle velocity gradients interpolated from the grid.
+    [[nodiscard]] ConstArrayView1<MatrixType> VelocityGradients() const;
+
+    //! Returns particle velocity gradients interpolated from the grid.
+    [[nodiscard]] ArrayView1<MatrixType> VelocityGradients();
+
+    //!
+    //! Transfers updated grid velocities and advances fluid kinematic state.
+    //!
+    //! The volume ratio is updated by the determinant of the incremental
+    //! deformation gradient, I + dt * velocityGradient.
+    //!
+    void TransferFromGridToParticles(double timeStepInSeconds);
+
+ private:
+    void ResetFluidState();
+
+    Array1<double> m_volumeRatios;
+    Array1<MatrixType> m_velocityGradients;
+};
+
+//! 2-D weakly compressible MPM fluid system data.
+using MPMFluidSystemData2 = MPMFluidSystemData<2>;
+
+//! 3-D weakly compressible MPM fluid system data.
+using MPMFluidSystemData3 = MPMFluidSystemData<3>;
+
+//! Shared pointer type of MPMFluidSystemData2.
+using MPMFluidSystemData2Ptr = std::shared_ptr<MPMFluidSystemData2>;
+
+//! Shared pointer type of MPMFluidSystemData3.
+using MPMFluidSystemData3Ptr = std::shared_ptr<MPMFluidSystemData3>;
+}  // namespace CubbyFlow
+
+#include <Core/Particle/MPM/MPMFluidSystemData-Impl.hpp>
+
+#endif
diff --git a/Tests/UnitTests/MPMFluidSystemDataTests.cpp b/Tests/UnitTests/MPMFluidSystemDataTests.cpp
new file mode 100644
index 0000000000..3adb137a59
--- /dev/null
+++ b/Tests/UnitTests/MPMFluidSystemDataTests.cpp
@@ -0,0 +1,474 @@
+// 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.
+
+#include "gtest/gtest.h"
+
+#include <Core/Particle/MPM/MPMFluidSystemData.hpp>
+
+#include <limits>
+#include <type_traits>
+
+using namespace CubbyFlow;
+
+namespace
+{
+template <typename T>
+concept HasDeformationStates = requires(T& data) { data.DeformationStates(); };
+
+template <typename T>
+concept HasNoArgumentGridToParticleTransfer =
+    requires(T& data) { data.TransferFromGridToParticles(); };
+
+static_assert(!HasDeformationStates<MPMFluidSystemData2>);
+static_assert(!HasDeformationStates<MPMFluidSystemData3>);
+static_assert(!HasNoArgumentGridToParticleTransfer<MPMTransferSystemData<2>>);
+static_assert(!HasNoArgumentGridToParticleTransfer<MPMTransferSystemData<3>>);
+static_assert(HasNoArgumentGridToParticleTransfer<MPMSystemData2>);
+static_assert(HasNoArgumentGridToParticleTransfer<MPMSystemData3>);
+static_assert(std::is_same_v<MPMSystemData2::Base, ParticleSystemData2>);
+static_assert(std::is_same_v<MPMSystemData3::Base, ParticleSystemData3>);
+
+template <size_t N>
+void ExpectParticleStateResizes()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+
+    EXPECT_DOUBLE_EQ(data.VolumeRatios()[0], 1.0);
+    EXPECT_TRUE(data.VelocityGradients()[0].IsSimilar(Matrix<double, N, N>{}));
+
+    data.VolumeRatios()[0] = 2.0;
+    data.VelocityGradients()[0](0, 0) = 3.0;
+    data.AddParticle(Vector<double, N>::MakeConstant(0.25));
+
+    ASSERT_EQ(data.VolumeRatios().Length(), 2u);
+    ASSERT_EQ(data.VelocityGradients().Length(), 2u);
+    EXPECT_DOUBLE_EQ(data.VolumeRatios()[0], 2.0);
+    EXPECT_DOUBLE_EQ(data.VolumeRatios()[1], 1.0);
+    EXPECT_DOUBLE_EQ(data.VelocityGradients()[0](0, 0), 3.0);
+    EXPECT_TRUE(data.VelocityGradients()[1].IsSimilar(Matrix<double, N, N>{}));
+}
+
+template <size_t N>
+void DirtyFluidState(MPMFluidSystemData<N>* data)
+{
+    data->VolumeRatios()[0] = 2.0;
+    data->VelocityGradients()[0](0, 0) = 3.0;
+}
+
+template <size_t N>
+void ExpectFluidStateReset(const MPMFluidSystemData<N>& data)
+{
+    ASSERT_EQ(data.VolumeRatios().Length(), 3u);
+    ASSERT_EQ(data.VelocityGradients().Length(), 3u);
+
+    for (size_t i = 0; i < 3; ++i)
+    {
+        EXPECT_DOUBLE_EQ(data.VolumeRatios()[i], 1.0);
+        EXPECT_TRUE(
+            data.VelocityGradients()[i].IsSimilar(Matrix<double, N, N>{}));
+    }
+}
+
+template <size_t N>
+void ExpectBaseSetResetsFluidState()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(2),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    ParticleSystemData<N> source{ 3 };
+
+    DirtyFluidState(&data);
+    data.Set(source);
+    ExpectFluidStateReset(data);
+}
+
+template <size_t N>
+void ExpectBaseDeserializeResetsFluidState()
+{
+    ParticleSystemData<N> source{ 3 };
+
+    std::vector<uint8_t> buffer;
+    source.Serialize(&buffer);
+
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(2),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+
+    DirtyFluidState(&data);
+    data.Deserialize(buffer);
+    ExpectFluidStateReset(data);
+}
+
+template <size_t N>
+void ExpectParticleToGridConservation()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                2 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.25);
+    data.Positions()[1] = Vector<double, N>::MakeConstant(2.25);
+    data.Velocities()[0] = Vector<double, N>::MakeConstant(2.0);
+    data.Velocities()[1] = Vector<double, N>::MakeConstant(-1.0);
+    data.ParticleMasses()[0] = 2.0;
+    data.ParticleMasses()[1] = 3.0;
+
+    data.TransferFromParticlesToGrid();
+
+    double totalMass = 0.0;
+    Vector<double, N> totalMomentum;
+    data.GridMass().ForEachDataPointIndex(
+        [&data, &totalMass, &totalMomentum](const Vector<size_t, N>& index) {
+            const double mass = data.GridMass()(index);
+            totalMass += mass;
+            totalMomentum += mass * data.GridVelocities()(index);
+        });
+
+    EXPECT_NEAR(totalMass, 5.0, 1e-11);
+    EXPECT_TRUE(
+        totalMomentum.IsSimilar(Vector<double, N>::MakeConstant(1.0), 1e-11));
+}
+
+template <size_t N>
+void ExpectParticleToGridFailureDoesNotCommit()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                2 };
+    const Vector<double, N> gridVelocity = Vector<double, N>::MakeConstant(3.0);
+    const Vector<double, N> oldGridVelocity =
+        Vector<double, N>::MakeConstant(4.0);
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.25);
+    data.Positions()[1] =
+        Vector<double, N>::MakeConstant(std::numeric_limits<double>::max());
+    data.GridMass().Fill(2.0, ExecutionPolicy::Serial);
+    data.GridVelocities().Fill(gridVelocity, ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Fill(oldGridVelocity,
+                                           ExecutionPolicy::Serial);
+
+    EXPECT_THROW(data.TransferFromParticlesToGrid(), std::invalid_argument);
+
+    data.GridMass().ForEachDataPointIndex(
+        [&data, &gridVelocity,
+         &oldGridVelocity](const Vector<size_t, N>& index) {
+            EXPECT_DOUBLE_EQ(data.GridMass()(index), 2.0);
+            EXPECT_TRUE(data.GridVelocities()(index).IsSimilar(gridVelocity));
+            EXPECT_TRUE(data.GridVelocitiesBeforeUpdate()(index).IsSimilar(
+                oldGridVelocity));
+        });
+}
+
+template <size_t N>
+void ExpectOverflowingParticleToGridTransferDoesNotCommit()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    const Vector<double, N> gridVelocity = Vector<double, N>::MakeConstant(3.0);
+    const Vector<double, N> oldGridVelocity =
+        Vector<double, N>::MakeConstant(4.0);
+
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.25);
+    data.Velocities()[0] =
+        Vector<double, N>::MakeConstant(std::numeric_limits<double>::max());
+    data.ParticleMasses()[0] = std::numeric_limits<double>::max();
+    data.GridMass().Fill(2.0, ExecutionPolicy::Serial);
+    data.GridVelocities().Fill(gridVelocity, ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Fill(oldGridVelocity,
+                                           ExecutionPolicy::Serial);
+
+    EXPECT_THROW(data.TransferFromParticlesToGrid(), std::invalid_argument);
+
+    data.GridMass().ForEachDataPointIndex(
+        [&data, &gridVelocity,
+         &oldGridVelocity](const Vector<size_t, N>& index) {
+            EXPECT_DOUBLE_EQ(data.GridMass()(index), 2.0);
+            EXPECT_TRUE(data.GridVelocities()(index).IsSimilar(gridVelocity));
+            EXPECT_TRUE(data.GridVelocitiesBeforeUpdate()(index).IsSimilar(
+                oldGridVelocity));
+        });
+}
+
+template <size_t N>
+void ExpectUniformVelocityRoundTrip()
+{
+    const Vector<double, N> velocity = Vector<double, N>::MakeConstant(2.5);
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(6),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                2 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(2.25);
+    data.Positions()[1] = Vector<double, N>::MakeConstant(3.25);
+    data.Velocities()[0] = velocity;
+    data.Velocities()[1] = velocity;
+
+    data.TransferFromParticlesToGrid();
+    data.TransferFromGridToParticles(0.1);
+
+    EXPECT_TRUE(data.Velocities()[0].IsSimilar(velocity, 1e-12));
+    EXPECT_TRUE(data.Velocities()[1].IsSimilar(velocity, 1e-12));
+    EXPECT_NEAR(data.VelocityGradients()[0].AbsMax(), 0.0, 1e-12);
+    EXPECT_NEAR(data.VelocityGradients()[1].AbsMax(), 0.0, 1e-12);
+    EXPECT_NEAR(data.VolumeRatios()[0], 1.0, 1e-12);
+    EXPECT_NEAR(data.VolumeRatios()[1], 1.0, 1e-12);
+}
+
+template <size_t N>
+void ExpectUniformDilation(double rate, double expectedVolumeRatio)
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(6),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(2.25);
+    data.GridVelocities().Fill(
+        [rate](const Vector<double, N>& position) { return rate * position; },
+        ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Set(data.GridVelocities());
+
+    data.TransferFromGridToParticles(1.0);
+
+    const auto& gradient = data.VelocityGradients()[0];
+
+    for (size_t row = 0; row < N; ++row)
+    {
+        for (size_t column = 0; column < N; ++column)
+        {
+            EXPECT_NEAR(gradient(row, column), row == column ? rate : 0.0,
+                        1e-12);
+        }
+    }
+
+    EXPECT_NEAR(data.VolumeRatios()[0], expectedVolumeRatio, 1e-12);
+}
+
+template <size_t N>
+void FillCrossVelocityField(MPMFluidSystemData<N>* data)
+{
+    data->GridVelocities().Fill(
+        [](const Vector<double, N>& position) {
+            Vector<double, N> velocity;
+            velocity[0] = 2.0 * position[1];
+            velocity[1] = 3.0 * position[0];
+            return velocity;
+        },
+        ExecutionPolicy::Serial);
+    data->GridVelocitiesBeforeUpdate().Set(data->GridVelocities());
+}
+
+template <size_t N>
+void ExpectVelocityGradientPreservesComponentOrientation()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(6),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(2.25);
+
+    FillCrossVelocityField(&data);
+    data.TransferFromGridToParticles(0.0);
+
+    const auto& gradient = data.VelocityGradients()[0];
+
+    for (size_t row = 0; row < N; ++row)
+    {
+        for (size_t column = 0; column < N; ++column)
+        {
+            double expected = 0.0;
+
+            if (row == 0 && column == 1)
+            {
+                expected = 2.0;
+            }
+            else if (row == 1 && column == 0)
+            {
+                expected = 3.0;
+            }
+
+            EXPECT_NEAR(gradient(row, column), expected, 1e-12);
+        }
+    }
+}
+
+template <size_t N>
+void ExpectVolumeRatioMultipliesNonDiagonalIncrement()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(6),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(2.25);
+    data.VolumeRatios()[0] = 2.0;
+
+    FillCrossVelocityField(&data);
+
+    data.TransferFromGridToParticles(0.1);
+    EXPECT_NEAR(data.VolumeRatios()[0], 1.88, 1e-12);
+
+    data.TransferFromGridToParticles(0.1);
+    EXPECT_NEAR(data.VolumeRatios()[0], 1.7672, 1e-12);
+}
+
+template <size_t N>
+void ExpectEmptyAndBoundaryTransfersAreSafe()
+{
+    MPMFluidSystemData<N> empty;
+    EXPECT_NO_THROW(empty.TransferFromParticlesToGrid());
+    EXPECT_NO_THROW(empty.TransferFromGridToParticles(0.1));
+
+    MPMFluidSystemData<N> boundary{ Vector<size_t, N>::MakeConstant(2),
+                                    Vector<double, N>::MakeConstant(1.0),
+                                    {},
+                                    1 };
+    boundary.Velocities()[0] = Vector<double, N>::MakeConstant(1.0);
+
+    EXPECT_NO_THROW(boundary.TransferFromParticlesToGrid());
+    EXPECT_NO_THROW(boundary.TransferFromGridToParticles(0.1));
+    EXPECT_TRUE(boundary.Velocities()[0].IsSimilar(
+        Vector<double, N>::MakeConstant(1.0), 1e-12));
+    EXPECT_NEAR(boundary.VolumeRatios()[0], 1.0, 1e-12);
+}
+
+template <size_t N>
+void ExpectInvalidFluidStateIsNotCommitted()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    const Vector<double, N> velocity = Vector<double, N>::MakeConstant(4.0);
+
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.25);
+    data.Velocities()[0] = velocity;
+    data.VolumeRatios()[0] = 0.0;
+    data.GridVelocities().Fill(Vector<double, N>::MakeConstant(2.0),
+                               ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Fill(Vector<double, N>::MakeConstant(1.0),
+                                           ExecutionPolicy::Serial);
+
+    EXPECT_THROW(data.TransferFromGridToParticles(0.1), std::invalid_argument);
+    EXPECT_TRUE(data.Velocities()[0].IsSimilar(velocity));
+    EXPECT_DOUBLE_EQ(data.VolumeRatios()[0], 0.0);
+    EXPECT_TRUE(data.VelocityGradients()[0].IsSimilar(Matrix<double, N, N>{}));
+    EXPECT_THROW(data.TransferFromGridToParticles(-0.1), std::invalid_argument);
+}
+
+template <size_t N>
+void ExpectOverflowingGridTransferIsNotCommitted()
+{
+    MPMFluidSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                                Vector<double, N>::MakeConstant(1.0),
+                                {},
+                                1 };
+    const Vector<double, N> velocity = Vector<double, N>::MakeConstant(4.0);
+    const double large = 0.6 * std::numeric_limits<double>::max();
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.25);
+    data.Velocities()[0] = velocity;
+    data.VolumeRatios()[0] = 2.0;
+    data.VelocityGradients()[0](0, 0) = 3.0;
+    data.GridVelocities().Fill(Vector<double, N>::MakeConstant(large),
+                               ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Fill(
+        Vector<double, N>::MakeConstant(-large), ExecutionPolicy::Serial);
+    data.SetFLIPBlendingFactor(1.0);
+
+    EXPECT_THROW(data.TransferFromGridToParticles(0.0), std::invalid_argument);
+    EXPECT_TRUE(data.Velocities()[0].IsSimilar(velocity));
+    EXPECT_DOUBLE_EQ(data.VolumeRatios()[0], 2.0);
+    EXPECT_DOUBLE_EQ(data.VelocityGradients()[0](0, 0), 3.0);
+}
+}  // namespace
+
+TEST(MPMFluidSystemData, ParticleStateResizes)
+{
+    ExpectParticleStateResizes<2>();
+    ExpectParticleStateResizes<3>();
+}
+
+TEST(MPMFluidSystemData, ParticleToGridConservesMassAndMomentum)
+{
+    ExpectParticleToGridConservation<2>();
+    ExpectParticleToGridConservation<3>();
+}
+
+TEST(MPMFluidSystemData, ParticleToGridFailureDoesNotCommit)
+{
+    ExpectParticleToGridFailureDoesNotCommit<2>();
+    ExpectParticleToGridFailureDoesNotCommit<3>();
+}
+
+TEST(MPMFluidSystemData, OverflowingParticleToGridTransferDoesNotCommit)
+{
+    ExpectOverflowingParticleToGridTransferDoesNotCommit<2>();
+    ExpectOverflowingParticleToGridTransferDoesNotCommit<3>();
+}
+
+TEST(MPMFluidSystemData, BaseSetResetsFluidState)
+{
+    ExpectBaseSetResetsFluidState<2>();
+    ExpectBaseSetResetsFluidState<3>();
+}
+
+TEST(MPMFluidSystemData, BaseDeserializeResetsFluidState)
+{
+    ExpectBaseDeserializeResetsFluidState<2>();
+    ExpectBaseDeserializeResetsFluidState<3>();
+}
+
+TEST(MPMFluidSystemData, UniformVelocityRoundTrip)
+{
+    ExpectUniformVelocityRoundTrip<2>();
+    ExpectUniformVelocityRoundTrip<3>();
+}
+
+TEST(MPMFluidSystemData, UniformExpansionAndCompression)
+{
+    ExpectUniformDilation<2>(0.1, 1.21);
+    ExpectUniformDilation<3>(0.1, 1.331);
+    ExpectUniformDilation<2>(-0.1, 0.81);
+    ExpectUniformDilation<3>(-0.1, 0.729);
+}
+
+TEST(MPMFluidSystemData, VelocityGradientPreservesComponentOrientation)
+{
+    ExpectVelocityGradientPreservesComponentOrientation<2>();
+    ExpectVelocityGradientPreservesComponentOrientation<3>();
+}
+
+TEST(MPMFluidSystemData, VolumeRatioMultipliesNonDiagonalIncrement)
+{
+    ExpectVolumeRatioMultipliesNonDiagonalIncrement<2>();
+    ExpectVolumeRatioMultipliesNonDiagonalIncrement<3>();
+}
+
+TEST(MPMFluidSystemData, EmptyAndBoundaryTransfersAreSafe)
+{
+    ExpectEmptyAndBoundaryTransfersAreSafe<2>();
+    ExpectEmptyAndBoundaryTransfersAreSafe<3>();
+}
+
+TEST(MPMFluidSystemData, InvalidFluidStateIsNotCommitted)
+{
+    ExpectInvalidFluidStateIsNotCommitted<2>();
+    ExpectInvalidFluidStateIsNotCommitted<3>();
+}
+
+TEST(MPMFluidSystemData, OverflowingGridTransferIsNotCommitted)
+{
+    ExpectOverflowingGridTransferIsNotCommitted<2>();
+    ExpectOverflowingGridTransferIsNotCommitted<3>();
+}

From 258eb1c4b69a92c2ea38c4b76a8224df8bb72d7b Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Mon, 24 Aug 2026 21:06:29 +0900
Subject: [PATCH 3/8] feat: expose fluid MPM state to Python

---
 .../API/Python/Particle/MPM/MPMSystemData.hpp |   2 +
 .../API/Python/Particle/MPM/MPMSystemData.cpp | 125 +++++++++++++++---
 Sources/API/Python/main.cpp                   |   2 +
 Tests/PythonTests/test_mpm_system_data.py     |  55 ++++++++
 4 files changed, 165 insertions(+), 19 deletions(-)

diff --git a/Includes/API/Python/Particle/MPM/MPMSystemData.hpp b/Includes/API/Python/Particle/MPM/MPMSystemData.hpp
index 5fc8a52d17..214a95187e 100644
--- a/Includes/API/Python/Particle/MPM/MPMSystemData.hpp
+++ b/Includes/API/Python/Particle/MPM/MPMSystemData.hpp
@@ -15,5 +15,7 @@
 
 void AddMPMSystemData2(pybind11::module& m);
 void AddMPMSystemData3(pybind11::module& m);
+void AddMPMFluidSystemData2(pybind11::module& m);
+void AddMPMFluidSystemData3(pybind11::module& m);
 
 #endif
diff --git a/Sources/API/Python/Particle/MPM/MPMSystemData.cpp b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
index f50e63444d..dc34f3de63 100644
--- a/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
+++ b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
@@ -9,6 +9,7 @@
 // property of any third parties.
 
 #include <API/Python/Particle/MPM/MPMSystemData.hpp>
+#include <Core/Particle/MPM/MPMFluidSystemData.hpp>
 #include <Core/Particle/MPM/MPMSystemData.hpp>
 
 #include <pybind11/numpy.h>
@@ -23,6 +24,11 @@ template <size_t N>
 using MPMClass = py::class_<MPMSystemData<N>, std::shared_ptr<MPMSystemData<N>>,
                             ParticleSystemData<N>>;
 
+template <size_t N>
+using MPMFluidClass =
+    py::class_<MPMFluidSystemData<N>, std::shared_ptr<MPMFluidSystemData<N>>,
+               ParticleSystemData<N>>;
+
 template <size_t N>
 py::array_t<double> MatrixView(SnowDeformationState<N>& state,
                                Matrix<double, N, N>& matrix)
@@ -48,15 +54,38 @@ void AddSnowDeformationState(py::module& m, const char* name)
         });
 }
 
-template <size_t N>
-py::array_t<double> ScalarView(MPMSystemData<N>& instance,
-                               ArrayView1<double> view)
+template <typename Data>
+py::array_t<double> ScalarView(Data& instance, ArrayView1<double> view)
 {
     return py::array_t<double>({ static_cast<py::ssize_t>(view.Length()) },
                                { static_cast<py::ssize_t>(sizeof(double)) },
                                view.data(), py::cast(&instance));
 }
 
+template <size_t N>
+py::array_t<double> FluidMatrixView(MPMFluidSystemData<N>& instance,
+                                    Matrix<double, N, N>& matrix)
+{
+    return py::array_t<double>(
+        { static_cast<py::ssize_t>(N), static_cast<py::ssize_t>(N) },
+        { static_cast<py::ssize_t>(sizeof(double) * N),
+          static_cast<py::ssize_t>(sizeof(double)) },
+        matrix.data(), py::cast(&instance));
+}
+
+template <size_t N>
+py::list VelocityGradients(MPMFluidSystemData<N>& instance)
+{
+    py::list result;
+
+    for (auto& gradient : instance.VelocityGradients())
+    {
+        result.append(FluidMatrixView(instance, gradient));
+    }
+
+    return result;
+}
+
 template <size_t N>
 py::list DeformationStates(MPMSystemData<N>& instance)
 {
@@ -87,34 +116,57 @@ void BindParticleState(MPMClass<N>& cls)
 }
 
 template <size_t N>
-void BindGridState(MPMClass<N>& cls)
+void BindParticleState(MPMFluidClass<N>& cls)
+{
+    cls.def_property_readonly("particleMasses",
+                              [](MPMFluidSystemData<N>& instance) {
+                                  return ScalarView(instance,
+                                                    instance.ParticleMasses());
+                              })
+        .def_property_readonly("initialVolumes",
+                               [](MPMFluidSystemData<N>& instance) {
+                                   return ScalarView(instance,
+                                                     instance.InitialVolumes());
+                               })
+        .def_property_readonly("volumeRatios",
+                               [](MPMFluidSystemData<N>& instance) {
+                                   return ScalarView(instance,
+                                                     instance.VolumeRatios());
+                               })
+        .def_property_readonly("velocityGradients", &VelocityGradients<N>);
+}
+
+template <size_t N, typename DataClass, typename Data>
+void BindGridState(DataClass& cls)
 {
-    cls.def("ResizeGrid", &MPMSystemData<N>::ResizeGrid)
+    cls.def("ResizeGrid",
+            [](Data& instance, const Vector<size_t, N>& resolution,
+               const Vector<double, N>& gridSpacing,
+               const Vector<double, N>& gridOrigin) {
+                instance.ResizeGrid(resolution, gridSpacing, gridOrigin);
+            })
         .def_property_readonly(
             "gridMass",
-            [](MPMSystemData<N>& instance) -> auto& {
-                return instance.GridMass();
-            },
+            [](Data& instance) -> auto& { return instance.GridMass(); },
             py::return_value_policy::reference_internal)
         .def_property_readonly(
             "gridVelocities",
-            [](MPMSystemData<N>& instance) -> auto& {
-                return instance.GridVelocities();
-            },
+            [](Data& instance) -> auto& { return instance.GridVelocities(); },
             py::return_value_policy::reference_internal)
         .def_property_readonly(
             "gridVelocitiesBeforeUpdate",
-            [](MPMSystemData<N>& instance) -> auto& {
+            [](Data& instance) -> auto& {
                 return instance.GridVelocitiesBeforeUpdate();
             },
             py::return_value_policy::reference_internal)
-        .def_property("flipBlendingFactor",
-                      &MPMSystemData<N>::FLIPBlendingFactor,
-                      &MPMSystemData<N>::SetFLIPBlendingFactor)
+        .def_property(
+            "flipBlendingFactor",
+            [](const Data& instance) { return instance.FLIPBlendingFactor(); },
+            [](Data& instance, double factor) {
+                instance.SetFLIPBlendingFactor(factor);
+            })
         .def("TransferFromParticlesToGrid",
-             &MPMSystemData<N>::TransferFromParticlesToGrid)
-        .def("TransferFromGridToParticles",
-             &MPMSystemData<N>::TransferFromGridToParticles);
+             [](Data& instance) { instance.TransferFromParticlesToGrid(); });
 }
 
 template <size_t N>
@@ -122,6 +174,7 @@ void AddMPMSystemData(py::module& m, const char* className,
                       const char* stateName)
 {
     AddSnowDeformationState<N>(m, stateName);
+
     MPMClass<N> cls(m, className);
     cls.def(py::init<const Vector<size_t, N>&, const Vector<double, N>&,
                      const Vector<double, N>&, size_t>(),
@@ -129,8 +182,32 @@ void AddMPMSystemData(py::module& m, const char* className,
             py::arg("gridSpacing") = Vector<double, N>::MakeConstant(1.0),
             py::arg("gridOrigin") = Vector<double, N>{},
             py::arg("numberOfParticles") = 0);
+
     BindParticleState(cls);
-    BindGridState(cls);
+    BindGridState<N, MPMClass<N>, MPMSystemData<N>>(cls);
+
+    cls.def("TransferFromGridToParticles", [](MPMSystemData<N>& instance) {
+        instance.TransferFromGridToParticles();
+    });
+}
+
+template <size_t N>
+void AddMPMFluidSystemData(py::module& m, const char* className)
+{
+    MPMFluidClass<N> cls(m, className);
+    cls.def(py::init<const Vector<size_t, N>&, const Vector<double, N>&,
+                     const Vector<double, N>&, size_t>(),
+            py::arg("resolution") = Vector<size_t, N>::MakeConstant(1),
+            py::arg("gridSpacing") = Vector<double, N>::MakeConstant(1.0),
+            py::arg("gridOrigin") = Vector<double, N>{},
+            py::arg("numberOfParticles") = 0);
+
+    BindParticleState(cls);
+    BindGridState<N, MPMFluidClass<N>, MPMFluidSystemData<N>>(cls);
+
+    cls.def("TransferFromGridToParticles",
+            &MPMFluidSystemData<N>::TransferFromGridToParticles,
+            py::arg("timeStepInSeconds"));
 }
 }  // namespace
 
@@ -143,3 +220,13 @@ void AddMPMSystemData3(pybind11::module& m)
 {
     AddMPMSystemData<3>(m, "MPMSystemData3", "SnowDeformationState3");
 }
+
+void AddMPMFluidSystemData2(pybind11::module& m)
+{
+    AddMPMFluidSystemData<2>(m, "MPMFluidSystemData2");
+}
+
+void AddMPMFluidSystemData3(pybind11::module& m)
+{
+    AddMPMFluidSystemData<3>(m, "MPMFluidSystemData3");
+}
diff --git a/Sources/API/Python/main.cpp b/Sources/API/Python/main.cpp
index 27f6303537..100d9c12b4 100644
--- a/Sources/API/Python/main.cpp
+++ b/Sources/API/Python/main.cpp
@@ -222,6 +222,8 @@ PYBIND11_MODULE(pyCubbyFlow, m)
     AddParticleSystemData3(m);
     AddMPMSystemData2(m);
     AddMPMSystemData3(m);
+    AddMPMFluidSystemData2(m);
+    AddMPMFluidSystemData3(m);
     AddSPHSystemData2(m);
     AddSPHSystemData3(m);
 
diff --git a/Tests/PythonTests/test_mpm_system_data.py b/Tests/PythonTests/test_mpm_system_data.py
index f6c391fff7..50e58c690a 100644
--- a/Tests/PythonTests/test_mpm_system_data.py
+++ b/Tests/PythonTests/test_mpm_system_data.py
@@ -63,3 +63,58 @@ def test_mpm_system_data_api(
         assert grid.resolution == resized_grid[0]
         assert grid.gridSpacing == resized_grid[1]
         assert grid.gridOrigin == resized_grid[2]
+
+
+@pytest.mark.parametrize(
+    "data_type,resolution,spacing,origin,grid_velocity,dimension",
+    [
+        (
+            pyCubbyFlow.MPMFluidSystemData2,
+            pyCubbyFlow.Vector2UZ(3, 3),
+            pyCubbyFlow.Vector2D(1.0, 1.0),
+            pyCubbyFlow.Vector2D(),
+            (1.0, 1.0),
+            2,
+        ),
+        (
+            pyCubbyFlow.MPMFluidSystemData3,
+            pyCubbyFlow.Vector3UZ(3, 3, 3),
+            pyCubbyFlow.Vector3D(1.0, 1.0, 1.0),
+            pyCubbyFlow.Vector3D(),
+            (1.0, 1.0, 1.0),
+            3,
+        ),
+    ],
+)
+def test_mpm_fluid_system_data_api(
+    data_type,
+    resolution,
+    spacing,
+    origin,
+    grid_velocity,
+    dimension,
+):
+    data = data_type(resolution, spacing, origin, 1)
+
+    ratios = np.asarray(data.volumeRatios)
+    gradients = data.velocityGradients
+    assert ratios.shape == (1,)
+    assert len(gradients) == 1
+    assert gradients[0].shape == (dimension, dimension)
+    assert not hasattr(data, "deformationStates")
+
+    ratios[0] = 2.0
+    gradients[0][0, 0] = 3.0
+    assert data.volumeRatios[0] == 2.0
+    assert data.velocityGradients[0][0, 0] == 3.0
+    ratios[0] = 1.0
+    gradients[0][0, 0] = 0.0
+
+    data.TransferFromParticlesToGrid()
+    data.gridVelocities.Fill(grid_velocity)
+    data.TransferFromGridToParticles(0.1)
+
+    assert data.volumeRatios[0] == pytest.approx(1.0)
+    assert np.asarray(data.velocityGradients[0]) == pytest.approx(
+        np.zeros((dimension, dimension))
+    )

From ec11886cdb0460c638f08173e08e30ed6cff72a9 Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Mon, 24 Aug 2026 22:51:22 +0900
Subject: [PATCH 4/8] fix: copy velocity gradient in tests to prevent dangling
 reference

---
 Tests/UnitTests/MPMFluidSystemDataTests.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Tests/UnitTests/MPMFluidSystemDataTests.cpp b/Tests/UnitTests/MPMFluidSystemDataTests.cpp
index 3adb137a59..365f40f4a2 100644
--- a/Tests/UnitTests/MPMFluidSystemDataTests.cpp
+++ b/Tests/UnitTests/MPMFluidSystemDataTests.cpp
@@ -242,7 +242,7 @@ void ExpectUniformDilation(double rate, double expectedVolumeRatio)
 
     data.TransferFromGridToParticles(1.0);
 
-    const auto& gradient = data.VelocityGradients()[0];
+    const auto gradient = data.VelocityGradients()[0];
 
     for (size_t row = 0; row < N; ++row)
     {
@@ -282,7 +282,7 @@ void ExpectVelocityGradientPreservesComponentOrientation()
     FillCrossVelocityField(&data);
     data.TransferFromGridToParticles(0.0);
 
-    const auto& gradient = data.VelocityGradients()[0];
+    const auto gradient = data.VelocityGradients()[0];
 
     for (size_t row = 0; row < N; ++row)
     {

From d7d8e64f31faebc457e27f377d973e3d39300612 Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Wed, 26 Aug 2026 11:02:54 +0900
Subject: [PATCH 5/8] fix: preserve MPM Python views during transfer

---
 .../Particle/MPM/MPMFluidSystemData-Impl.hpp  |  7 +-
 .../API/Python/Particle/MPM/MPMSystemData.cpp | 74 ++++++++++++-------
 Tests/PythonTests/test_mpm_system_data.py     | 54 ++++++++++++++
 3 files changed, 107 insertions(+), 28 deletions(-)

diff --git a/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
index f8787b746e..6305abea24 100644
--- a/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
+++ b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
@@ -11,6 +11,7 @@
 #ifndef CUBBYFLOW_MPM_FLUID_SYSTEM_DATA_IMPL_HPP
 #define CUBBYFLOW_MPM_FLUID_SYSTEM_DATA_IMPL_HPP
 
+#include <algorithm>
 #include <cmath>
 #include <stdexcept>
 
@@ -117,8 +118,10 @@ void MPMFluidSystemData<N>::TransferFromGridToParticles(
 
     this->TransferFromGridToParticlesUnchecked();
 
-    m_volumeRatios.Swap(nextVolumeRatios);
-    m_velocityGradients.Swap(nextVelocityGradients);
+    std::copy(nextVolumeRatios.begin(), nextVolumeRatios.end(),
+              m_volumeRatios.begin());
+    std::copy(nextVelocityGradients.begin(), nextVelocityGradients.end(),
+              m_velocityGradients.begin());
 }
 
 template <size_t N>
diff --git a/Sources/API/Python/Particle/MPM/MPMSystemData.cpp b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
index dc34f3de63..de8eb6217f 100644
--- a/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
+++ b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
@@ -20,6 +20,22 @@ namespace
 {
 namespace py = pybind11;
 
+constexpr char SCALAR_VIEW_DOC[] = R"pbdoc(
+Returns a mutable NumPy view into the particle data.
+
+The view is invalidated when an operation changes the particle count or
+replaces the particle state, such as Resize, AddParticle, AddParticles, Set,
+or Deserialize.
+)pbdoc";
+
+constexpr char MATRIX_VIEW_DOC[] = R"pbdoc(
+Returns mutable NumPy views into the per-particle matrices.
+
+The views are invalidated when an operation changes the particle count or
+replaces the particle state, such as Resize, AddParticle, AddParticles, Set,
+or Deserialize.
+)pbdoc";
+
 template <size_t N>
 using MPMClass = py::class_<MPMSystemData<N>, std::shared_ptr<MPMSystemData<N>>,
                             ParticleSystemData<N>>;
@@ -102,38 +118,44 @@ py::list DeformationStates(MPMSystemData<N>& instance)
 template <size_t N>
 void BindParticleState(MPMClass<N>& cls)
 {
-    cls.def_property_readonly("particleMasses",
-                              [](MPMSystemData<N>& instance) {
-                                  return ScalarView(instance,
-                                                    instance.ParticleMasses());
-                              })
-        .def_property_readonly("initialVolumes",
-                               [](MPMSystemData<N>& instance) {
-                                   return ScalarView(instance,
-                                                     instance.InitialVolumes());
-                               })
+    cls.def_property_readonly(
+           "particleMasses",
+           [](MPMSystemData<N>& instance) {
+               return ScalarView(instance, instance.ParticleMasses());
+           },
+           SCALAR_VIEW_DOC)
+        .def_property_readonly(
+            "initialVolumes",
+            [](MPMSystemData<N>& instance) {
+                return ScalarView(instance, instance.InitialVolumes());
+            },
+            SCALAR_VIEW_DOC)
         .def_property_readonly("deformationStates", &DeformationStates<N>);
 }
 
 template <size_t N>
 void BindParticleState(MPMFluidClass<N>& cls)
 {
-    cls.def_property_readonly("particleMasses",
-                              [](MPMFluidSystemData<N>& instance) {
-                                  return ScalarView(instance,
-                                                    instance.ParticleMasses());
-                              })
-        .def_property_readonly("initialVolumes",
-                               [](MPMFluidSystemData<N>& instance) {
-                                   return ScalarView(instance,
-                                                     instance.InitialVolumes());
-                               })
-        .def_property_readonly("volumeRatios",
-                               [](MPMFluidSystemData<N>& instance) {
-                                   return ScalarView(instance,
-                                                     instance.VolumeRatios());
-                               })
-        .def_property_readonly("velocityGradients", &VelocityGradients<N>);
+    cls.def_property_readonly(
+           "particleMasses",
+           [](MPMFluidSystemData<N>& instance) {
+               return ScalarView(instance, instance.ParticleMasses());
+           },
+           SCALAR_VIEW_DOC)
+        .def_property_readonly(
+            "initialVolumes",
+            [](MPMFluidSystemData<N>& instance) {
+                return ScalarView(instance, instance.InitialVolumes());
+            },
+            SCALAR_VIEW_DOC)
+        .def_property_readonly(
+            "volumeRatios",
+            [](MPMFluidSystemData<N>& instance) {
+                return ScalarView(instance, instance.VolumeRatios());
+            },
+            SCALAR_VIEW_DOC)
+        .def_property_readonly("velocityGradients", &VelocityGradients<N>,
+                               MATRIX_VIEW_DOC);
 }
 
 template <size_t N, typename DataClass, typename Data>
diff --git a/Tests/PythonTests/test_mpm_system_data.py b/Tests/PythonTests/test_mpm_system_data.py
index 50e58c690a..0663c5ca6b 100644
--- a/Tests/PythonTests/test_mpm_system_data.py
+++ b/Tests/PythonTests/test_mpm_system_data.py
@@ -118,3 +118,57 @@ def test_mpm_fluid_system_data_api(
     assert np.asarray(data.velocityGradients[0]) == pytest.approx(
         np.zeros((dimension, dimension))
     )
+
+
+@pytest.mark.parametrize(
+    "data_type,resolution,spacing,origin,position,velocity_field,dimension",
+    [
+        (
+            pyCubbyFlow.MPMFluidSystemData2,
+            pyCubbyFlow.Vector2UZ(6, 6),
+            pyCubbyFlow.Vector2D(1.0, 1.0),
+            pyCubbyFlow.Vector2D(),
+            (2.25, 2.25),
+            lambda pt: (0.2 * pt.x, 0.2 * pt.y),
+            2,
+        ),
+        (
+            pyCubbyFlow.MPMFluidSystemData3,
+            pyCubbyFlow.Vector3UZ(6, 6, 6),
+            pyCubbyFlow.Vector3D(1.0, 1.0, 1.0),
+            pyCubbyFlow.Vector3D(),
+            (2.25, 2.25, 2.25),
+            lambda pt: (0.2 * pt.x, 0.2 * pt.y, 0.2 * pt.z),
+            3,
+        ),
+    ],
+)
+def test_mpm_fluid_views_remain_valid_after_transfer(
+    data_type,
+    resolution,
+    spacing,
+    origin,
+    position,
+    velocity_field,
+    dimension,
+):
+    data = data_type(resolution, spacing, origin)
+    data.AddParticle(position)
+    data.gridVelocities.Fill(velocity_field)
+    data.gridVelocitiesBeforeUpdate.Fill(velocity_field)
+
+    ratios = np.asarray(data.volumeRatios)
+    gradients = data.velocityGradients
+    ratios[0] = 2.0
+
+    data.TransferFromGridToParticles(0.5)
+
+    assert ratios[0] == pytest.approx(2.0 * 1.1**dimension)
+    assert np.asarray(gradients[0]) == pytest.approx(
+        0.2 * np.identity(dimension)
+    )
+
+    ratios[0] = 3.0
+    gradients[0][0, 0] = 4.0
+    assert data.volumeRatios[0] == 3.0
+    assert data.velocityGradients[0][0, 0] == 4.0

From 4d15b1ba30bc73cd2c9314c2bb5d4b3d86fd6944 Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Wed, 26 Aug 2026 11:03:36 +0900
Subject: [PATCH 6/8] fix: avoid virtual resize in MPM constructor

---
 Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
index 3232811c66..4fb9ed4ecf 100644
--- a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
+++ b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
@@ -162,10 +162,11 @@ 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{}
+    : Base{ numberOfParticles },
+      m_particleMasses(numberOfParticles, Base::Mass()),
+      m_initialVolumes(numberOfParticles, 0.0)
 {
     ResizeGrid(resolution, gridSpacing, gridOrigin);
-    Resize(numberOfParticles);
 }
 
 template <size_t N>

From 38faa9d831172637fa27b7d5b2fa85ebda6bc1cf Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Wed, 26 Aug 2026 11:03:52 +0900
Subject: [PATCH 7/8] refactor: remove redundant MPM transfer wrapper

---
 Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp | 6 ------
 Includes/Core/Particle/MPM/MPMSystemData.hpp      | 5 +----
 2 files changed, 1 insertion(+), 10 deletions(-)

diff --git a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
index 4fb9ed4ecf..51be391c95 100644
--- a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
+++ b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
@@ -595,12 +595,6 @@ void MPMSystemData<N>::Set(const ParticleSystemData<N>& other)
     m_deformationStates.Fill(DeformationState{});
 }
 
-template <size_t N>
-void MPMSystemData<N>::TransferFromGridToParticles()
-{
-    TransferBase::TransferFromGridToParticles();
-}
-
 template <size_t N>
 ConstArrayView1<typename MPMSystemData<N>::DeformationState>
 MPMSystemData<N>::DeformationStates() const
diff --git a/Includes/Core/Particle/MPM/MPMSystemData.hpp b/Includes/Core/Particle/MPM/MPMSystemData.hpp
index 25b5c4407c..b4c2e85a2e 100644
--- a/Includes/Core/Particle/MPM/MPMSystemData.hpp
+++ b/Includes/Core/Particle/MPM/MPMSystemData.hpp
@@ -188,6 +188,7 @@ class MPMSystemData final : public MPMTransferSystemData<N>
     using Base = ParticleSystemData<N>;
     using TransferBase = MPMTransferSystemData<N>;
     using DeformationState = SnowDeformationState<N>;
+    using TransferBase::TransferFromGridToParticles;
 
     //! Constructs snow MPM state with a vertex-centered background grid.
     explicit MPMSystemData(
@@ -207,10 +208,6 @@ class MPMSystemData final : public MPMTransferSystemData<N>
     //! Copies inherited particle state and resets deformation state.
     void Set(const ParticleSystemData<N>& other) override;
 
-    //! Transfers updated background-grid velocities to particles.
-    //! Stencil nodes outside the finite grid are clamped to its boundary.
-    void TransferFromGridToParticles();
-
     //! Returns per-particle deformation states.
     [[nodiscard]] ConstArrayView1<DeformationState> DeformationStates() const;
 

From 242fe29519924566d7545f5bc291af54df96a06b Mon Sep 17 00:00:00 2001
From: Chris Ohk <utilforever@gmail.com>
Date: Wed, 26 Aug 2026 12:31:37 +0900
Subject: [PATCH 8/8] refactor: use std::ranges::copy for particle data
 transfer

---
 Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
index 6305abea24..d0bc66b08d 100644
--- a/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
+++ b/Includes/Core/Particle/MPM/MPMFluidSystemData-Impl.hpp
@@ -118,10 +118,8 @@ void MPMFluidSystemData<N>::TransferFromGridToParticles(
 
     this->TransferFromGridToParticlesUnchecked();
 
-    std::copy(nextVolumeRatios.begin(), nextVolumeRatios.end(),
-              m_volumeRatios.begin());
-    std::copy(nextVelocityGradients.begin(), nextVelocityGradients.end(),
-              m_velocityGradients.begin());
+    std::ranges::copy(nextVolumeRatios, m_volumeRatios.begin());
+    std::ranges::copy(nextVelocityGradients, m_velocityGradients.begin());
 }
 
 template <size_t N>