diff --git a/.github/workflows/ubuntu-sonarcloud.yml b/.github/workflows/ubuntu-sonarcloud.yml
index 043b71457b..7721e17e6a 100644
--- a/.github/workflows/ubuntu-sonarcloud.yml
+++ b/.github/workflows/ubuntu-sonarcloud.yml
@@ -51,7 +51,7 @@ jobs:
       - name: Run Unit Test
         run: |
           build/bin/UnitTests
-          gcovr --gcov-executable "gcov-${{ matrix.version }}" --root . --filter '^Sources/' --sonarqube build/coverage.xml build
+          gcovr --gcov-executable "gcov-${{ matrix.version }}" --root . --filter '^(Includes|Sources)/' --sonarqube build/coverage.xml build
       - name: SonarCloud Scan
         run: sonar-scanner -Dsonar.host.url=https://sonarcloud.io -Dsonar.organization=utilforever-github
         env:
diff --git a/Includes/API/Python/Particle/MPM/MPMSystemData.hpp b/Includes/API/Python/Particle/MPM/MPMSystemData.hpp
new file mode 100644
index 0000000000..5fc8a52d17
--- /dev/null
+++ b/Includes/API/Python/Particle/MPM/MPMSystemData.hpp
@@ -0,0 +1,19 @@
+// 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_PYTHON_MPM_SYSTEM_DATA_HPP
+#define CUBBYFLOW_PYTHON_MPM_SYSTEM_DATA_HPP
+
+#include <pybind11/pybind11.h>
+
+void AddMPMSystemData2(pybind11::module& m);
+void AddMPMSystemData3(pybind11::module& m);
+
+#endif
diff --git a/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
new file mode 100644
index 0000000000..0577661083
--- /dev/null
+++ b/Includes/Core/Particle/MPM/MPMSystemData-Impl.hpp
@@ -0,0 +1,497 @@
+// 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>
+
+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>
+MPMSystemData<N>::MPMSystemData(const Vector<size_t, N>& resolution,
+                                const Vector<double, N>& gridSpacing,
+                                const Vector<double, N>& gridOrigin,
+                                size_t numberOfParticles)
+    : Base{}
+{
+    ResizeGrid(resolution, gridSpacing, gridOrigin);
+    Resize(numberOfParticles);
+}
+
+template <size_t N>
+void MPMSystemData<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)
+{
+    Base::Deserialize(buffer);
+    ResetMPMState();
+}
+
+template <size_t N>
+void MPMSystemData<N>::Set(const ParticleSystemData<N>& other)
+{
+    Base::Set(other);
+    ResetMPMState();
+}
+
+template <size_t N>
+void MPMSystemData<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> MPMSystemData<N>::ParticleMasses() const
+{
+    return m_particleMasses.View();
+}
+
+template <size_t N>
+ArrayView1<double> MPMSystemData<N>::ParticleMasses()
+{
+    return m_particleMasses.View();
+}
+
+template <size_t N>
+ConstArrayView1<double> MPMSystemData<N>::InitialVolumes() const
+{
+    return m_initialVolumes.View();
+}
+
+template <size_t N>
+ArrayView1<double> MPMSystemData<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
+{
+    return m_gridMass;
+}
+
+template <size_t N>
+VertexCenteredScalarGrid<N>& MPMSystemData<N>::GridMass()
+{
+    return m_gridMass;
+}
+
+template <size_t N>
+const VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocities() const
+{
+    return m_gridVelocities;
+}
+
+template <size_t N>
+VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocities()
+{
+    return m_gridVelocities;
+}
+
+template <size_t N>
+const VertexCenteredVectorGrid<N>&
+MPMSystemData<N>::GridVelocitiesBeforeUpdate() const
+{
+    return m_gridVelocitiesBeforeUpdate;
+}
+
+template <size_t N>
+VertexCenteredVectorGrid<N>& MPMSystemData<N>::GridVelocitiesBeforeUpdate()
+{
+    return m_gridVelocitiesBeforeUpdate;
+}
+
+template <size_t N>
+double MPMSystemData<N>::FLIPBlendingFactor() const
+{
+    return m_flipBlendingFactor;
+}
+
+template <size_t N>
+void MPMSystemData<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 MPMSystemData<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.");
+        }
+    }
+
+    m_gridMass.Fill(0.0, ExecutionPolicy::Serial);
+    m_gridVelocities.Fill(Vector<double, N>{}, ExecutionPolicy::Serial);
+
+    const auto dataSize = m_gridMass.DataSize();
+    const auto gridSpacing = m_gridMass.GridSpacing();
+    const auto dataOrigin = m_gridMass.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];
+
+            m_gridMass(index) += mass;
+            m_gridVelocities(index) += mass * velocities[i];
+        }
+    }
+
+    m_gridMass.ForEachDataPointIndex([this](const Vector<size_t, N>& index) {
+        const double mass = m_gridMass(index);
+        if (mass > 0.0)
+        {
+            m_gridVelocities(index) /= mass;
+        }
+    });
+
+    m_gridVelocitiesBeforeUpdate.Set(m_gridVelocities);
+}
+
+template <size_t N>
+void MPMSystemData<N>::TransferFromGridToParticles()
+{
+    ValidateGridState();
+
+    const auto positions = this->Positions();
+    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.");
+            }
+        });
+
+    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;
+
+        velocities[i] = result;
+    }
+}
+
+template <size_t N>
+Vector<size_t, N> MPMSystemData<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 MPMSystemData<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 MPMSystemData<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 MPMSystemData<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 MPMSystemData<N>::ResetMPMState()
+{
+    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);
+}
+}  // namespace CubbyFlow
+
+#endif
diff --git a/Includes/Core/Particle/MPM/MPMSystemData.hpp b/Includes/Core/Particle/MPM/MPMSystemData.hpp
new file mode 100644
index 0000000000..008c6ca6bf
--- /dev/null
+++ b/Includes/Core/Particle/MPM/MPMSystemData.hpp
@@ -0,0 +1,191 @@
+// 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_HPP
+#define CUBBYFLOW_MPM_SYSTEM_DATA_HPP
+
+#include <Core/Grid/VertexCenteredScalarGrid.hpp>
+#include <Core/Grid/VertexCenteredVectorGrid.hpp>
+#include <Core/Particle/MPM/SnowConstitutiveModel.hpp>
+#include <Core/Particle/ParticleSystemData.hpp>
+
+#include <array>
+
+namespace CubbyFlow
+{
+//!
+//! \brief Tensor-product cubic B-spline interpolation kernel.
+//!
+template <size_t N>
+class CubicBSplineKernel final
+{
+ public:
+    static_assert(N == 2 || N == 3, "MPM supports only 2-D and 3-D.");
+
+    static constexpr size_t STENCIL_SIZE = N == 2 ? 16 : 64;
+
+    struct Entry
+    {
+        Vector<ssize_t, N> index;
+        double weight = 0.0;
+        Vector<double, N> gradient;
+    };
+
+    using Stencil = std::array<Entry, STENCIL_SIZE>;
+
+    [[nodiscard]] static double Weight(double x);
+
+    [[nodiscard]] static double Gradient(double x);
+
+    [[nodiscard]] static Stencil GetStencil(
+        const Vector<double, N>& position, const Vector<double, N>& gridSpacing,
+        const Vector<double, N>& dataOrigin);
+
+ private:
+    static void 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);
+
+    [[nodiscard]] static Entry GetStencilEntry(
+        const Vector<size_t, N>& offset, const Vector<double, N>& normalized,
+        const Vector<ssize_t, N>& firstIndex,
+        const Vector<double, N>& gridSpacing);
+};
+
+//!
+//! \brief N-D 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 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(
+        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 MPM attributes.
+    void Resize(size_t newNumberOfParticles) override;
+
+    //! Deserializes inherited particle state and resets MPM-specific state.
+    void Deserialize(const std::vector<uint8_t>& buffer) override;
+
+    //! Copies inherited particle state and resets MPM-specific state.
+    void Set(const ParticleSystemData<N>& other) override;
+
+    //! Resizes the background grid without changing particle state.
+    void ResizeGrid(const Vector<size_t, N>& resolution,
+                    const Vector<double, N>& gridSpacing,
+                    const Vector<double, N>& gridOrigin);
+
+    //! Returns per-particle masses.
+    [[nodiscard]] ConstArrayView1<double> ParticleMasses() const;
+
+    //! Returns per-particle masses.
+    [[nodiscard]] ArrayView1<double> ParticleMasses();
+
+    //! Returns per-particle initial volumes.
+    [[nodiscard]] ConstArrayView1<double> InitialVolumes() const;
+
+    //! 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;
+
+    //! Returns grid mass.
+    [[nodiscard]] VertexCenteredScalarGrid<N>& GridMass();
+
+    //! Returns current grid velocities.
+    [[nodiscard]] const VertexCenteredVectorGrid<N>& GridVelocities() const;
+
+    //! Returns current grid velocities.
+    [[nodiscard]] VertexCenteredVectorGrid<N>& GridVelocities();
+
+    //! Returns grid velocities before the grid update.
+    [[nodiscard]] const VertexCenteredVectorGrid<N>&
+    GridVelocitiesBeforeUpdate() const;
+
+    //! Returns grid velocities before the grid update.
+    [[nodiscard]] VertexCenteredVectorGrid<N>& GridVelocitiesBeforeUpdate();
+
+    //! Returns the FLIP fraction used for grid-to-particle transfer.
+    [[nodiscard]] double FLIPBlendingFactor() const;
+
+    //! Sets the FLIP fraction used for grid-to-particle transfer.
+    void SetFLIPBlendingFactor(double factor);
+
+    //! Transfers particle mass and momentum to the background grid.
+    //! Stencil nodes outside the finite grid are clamped to its boundary.
+    void TransferFromParticlesToGrid();
+
+    //! 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);
+
+    [[nodiscard]] static bool IsFinite(const Vector<double, N>& value);
+
+    static void ValidateGridParameters(const Vector<size_t, N>& resolution,
+                                       const Vector<double, N>& gridSpacing,
+                                       const Vector<double, N>& gridOrigin);
+
+    void ValidateGridState() const;
+
+    void ResetMPMState();
+
+    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;
+};
+
+//! 2-D material point method system data.
+using MPMSystemData2 = MPMSystemData<2>;
+
+//! 3-D material point method system data.
+using MPMSystemData3 = MPMSystemData<3>;
+
+//! Shared pointer type of MPMSystemData2.
+using MPMSystemData2Ptr = std::shared_ptr<MPMSystemData2>;
+
+//! Shared pointer type of MPMSystemData3.
+using MPMSystemData3Ptr = std::shared_ptr<MPMSystemData3>;
+}  // namespace CubbyFlow
+
+#include <Core/Particle/MPM/MPMSystemData-Impl.hpp>
+
+#endif
diff --git a/Includes/Core/Particle/ParticleSystemData.hpp b/Includes/Core/Particle/ParticleSystemData.hpp
index b4cf5a253a..8f05c2f674 100644
--- a/Includes/Core/Particle/ParticleSystemData.hpp
+++ b/Includes/Core/Particle/ParticleSystemData.hpp
@@ -86,7 +86,7 @@ class ParticleSystemData : public Serializable
     //!
     //! \param[in]  newNumberOfParticles    New number of particles.
     //!
-    void Resize(size_t newNumberOfParticles);
+    virtual void Resize(size_t newNumberOfParticles);
 
     //! Returns the number of particles.
     [[nodiscard]] size_t NumberOfParticles() const;
@@ -231,7 +231,7 @@ class ParticleSystemData : public Serializable
     void Deserialize(const std::vector<uint8_t>& buffer) override;
 
     //! Copies from other particle system data.
-    void Set(const ParticleSystemData& other);
+    virtual void Set(const ParticleSystemData& other);
 
  protected:
     template <size_t M = N>
@@ -284,4 +284,4 @@ using ParticleSystemData2Ptr = std::shared_ptr<ParticleSystemData2>;
 using ParticleSystemData3Ptr = std::shared_ptr<ParticleSystemData3>;
 }  // namespace CubbyFlow
 
-#endif
\ No newline at end of file
+#endif
diff --git a/Includes/Core/Particle/SPHSystemData.hpp b/Includes/Core/Particle/SPHSystemData.hpp
index 737f7772c1..1e7e103b17 100644
--- a/Includes/Core/Particle/SPHSystemData.hpp
+++ b/Includes/Core/Particle/SPHSystemData.hpp
@@ -35,6 +35,7 @@ class SPHSystemData : public ParticleSystemData<N>
     using Base::Positions;
     using Base::ScalarDataAt;
     using Base::Serialize;
+    using Base::Set;
 
     //! Constructs empty SPH system.
     SPHSystemData();
@@ -256,4 +257,4 @@ using SPHSystemData2Ptr = std::shared_ptr<SPHSystemData2>;
 using SPHSystemData3Ptr = std::shared_ptr<SPHSystemData3>;
 }  // namespace CubbyFlow
 
-#endif
\ No newline at end of file
+#endif
diff --git a/Sources/API/Python/Particle/MPM/MPMSystemData.cpp b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
new file mode 100644
index 0000000000..f50e63444d
--- /dev/null
+++ b/Sources/API/Python/Particle/MPM/MPMSystemData.cpp
@@ -0,0 +1,145 @@
+// 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 <API/Python/Particle/MPM/MPMSystemData.hpp>
+#include <Core/Particle/MPM/MPMSystemData.hpp>
+
+#include <pybind11/numpy.h>
+
+using namespace CubbyFlow;
+
+namespace
+{
+namespace py = pybind11;
+
+template <size_t N>
+using MPMClass = py::class_<MPMSystemData<N>, std::shared_ptr<MPMSystemData<N>>,
+                            ParticleSystemData<N>>;
+
+template <size_t N>
+py::array_t<double> MatrixView(SnowDeformationState<N>& state,
+                               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(&state));
+}
+
+template <size_t N>
+void AddSnowDeformationState(py::module& m, const char* name)
+{
+    using State = SnowDeformationState<N>;
+    py::class_<State>(m, name)
+        .def(py::init<>())
+        .def_property_readonly(
+            "elastic",
+            [](State& state) { return MatrixView(state, state.elastic); })
+        .def_property_readonly("plastic", [](State& state) {
+            return MatrixView(state, state.plastic);
+        });
+}
+
+template <size_t N>
+py::array_t<double> ScalarView(MPMSystemData<N>& 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::list DeformationStates(MPMSystemData<N>& instance)
+{
+    py::list result;
+    auto parent = py::cast(&instance, py::return_value_policy::reference);
+    for (auto& state : instance.DeformationStates())
+    {
+        result.append(
+            py::cast(&state, py::return_value_policy::reference, parent));
+    }
+    return result;
+}
+
+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());
+                               })
+        .def_property_readonly("deformationStates", &DeformationStates<N>);
+}
+
+template <size_t N>
+void BindGridState(MPMClass<N>& cls)
+{
+    cls.def("ResizeGrid", &MPMSystemData<N>::ResizeGrid)
+        .def_property_readonly(
+            "gridMass",
+            [](MPMSystemData<N>& instance) -> auto& {
+                return instance.GridMass();
+            },
+            py::return_value_policy::reference_internal)
+        .def_property_readonly(
+            "gridVelocities",
+            [](MPMSystemData<N>& instance) -> auto& {
+                return instance.GridVelocities();
+            },
+            py::return_value_policy::reference_internal)
+        .def_property_readonly(
+            "gridVelocitiesBeforeUpdate",
+            [](MPMSystemData<N>& instance) -> auto& {
+                return instance.GridVelocitiesBeforeUpdate();
+            },
+            py::return_value_policy::reference_internal)
+        .def_property("flipBlendingFactor",
+                      &MPMSystemData<N>::FLIPBlendingFactor,
+                      &MPMSystemData<N>::SetFLIPBlendingFactor)
+        .def("TransferFromParticlesToGrid",
+             &MPMSystemData<N>::TransferFromParticlesToGrid)
+        .def("TransferFromGridToParticles",
+             &MPMSystemData<N>::TransferFromGridToParticles);
+}
+
+template <size_t N>
+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>(),
+            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(cls);
+}
+}  // namespace
+
+void AddMPMSystemData2(pybind11::module& m)
+{
+    AddMPMSystemData<2>(m, "MPMSystemData2", "SnowDeformationState2");
+}
+
+void AddMPMSystemData3(pybind11::module& m)
+{
+    AddMPMSystemData<3>(m, "MPMSystemData3", "SnowDeformationState3");
+}
diff --git a/Sources/API/Python/Particle/SPH/SPHSystemData.cpp b/Sources/API/Python/Particle/SPH/SPHSystemData.cpp
index 8925eb78bd..6cbb7cfebc 100644
--- a/Sources/API/Python/Particle/SPH/SPHSystemData.cpp
+++ b/Sources/API/Python/Particle/SPH/SPHSystemData.cpp
@@ -89,7 +89,9 @@ void AddSPHSystemData2(pybind11::module& m)
              R"pbdoc(
 			Builds neighbor lists with kernel radius.
 		)pbdoc")
-        .def("Set", &SPHSystemData2::Set,
+        .def("Set",
+             static_cast<void (SPHSystemData2::*)(const SPHSystemData2&)>(
+                 &SPHSystemData2::Set),
              R"pbdoc(
 			Copies from other SPH system data.
 		)pbdoc");
@@ -169,8 +171,10 @@ void AddSPHSystemData3(pybind11::module& m)
              R"pbdoc(
 			Builds neighbor lists with kernel radius.
 		)pbdoc")
-        .def("Set", &SPHSystemData3::Set,
+        .def("Set",
+             static_cast<void (SPHSystemData3::*)(const SPHSystemData3&)>(
+                 &SPHSystemData3::Set),
              R"pbdoc(
 			Copies from other SPH system data.
 		)pbdoc");
-}
\ No newline at end of file
+}
diff --git a/Sources/API/Python/main.cpp b/Sources/API/Python/main.cpp
index b0e300d285..cce67f6902 100644
--- a/Sources/API/Python/main.cpp
+++ b/Sources/API/Python/main.cpp
@@ -54,6 +54,7 @@
 #include <API/Python/Grid/VertexCenteredScalarGrid.hpp>
 #include <API/Python/Grid/VertexCenteredVectorGrid.hpp>
 #include <API/Python/Math/Quaternion.hpp>
+#include <API/Python/Particle/MPM/MPMSystemData.hpp>
 #include <API/Python/Particle/ParticleSystemData.hpp>
 #include <API/Python/Particle/SPH/SPHSystemData.hpp>
 #include <API/Python/PointsToImplicit/AnisotropicPointsToImplicit.hpp>
@@ -218,6 +219,8 @@ PYBIND11_MODULE(pyCubbyFlow, m)
     AddGridSystemData3(m);
     AddParticleSystemData2(m);
     AddParticleSystemData3(m);
+    AddMPMSystemData2(m);
+    AddMPMSystemData3(m);
     AddSPHSystemData2(m);
     AddSPHSystemData3(m);
 
@@ -338,4 +341,4 @@ PYBIND11_MODULE(pyCubbyFlow, m)
 #else
     m.attr("__version__") = pybind11::str("Dev");
 #endif
-}
\ No newline at end of file
+}
diff --git a/Sources/Core/Particle/ParticleSystemData.cpp b/Sources/Core/Particle/ParticleSystemData.cpp
index b8a5aad118..1014895c6b 100644
--- a/Sources/Core/Particle/ParticleSystemData.cpp
+++ b/Sources/Core/Particle/ParticleSystemData.cpp
@@ -72,7 +72,7 @@ ParticleSystemData<N>::ParticleSystemData(size_t numberOfParticles)
         Vector<size_t, N>::MakeConstant(DEFAULT_HASH_GRID_RESOLUTION),
         2.0 * m_radius);
 
-    Resize(numberOfParticles);
+    ParticleSystemData<N>::Resize(numberOfParticles);
 }
 
 template <size_t N>
@@ -422,12 +422,14 @@ void ParticleSystemData<N>::Deserialize(const std::vector<uint8_t>& buffer)
         GetFlatbuffersParticleSystemData<N>::GetParticleSystemData(
             buffer.data());
     Deserialize(fbsParticleSystemData, *this);
+    Resize(NumberOfParticles());
 }
 
 template <size_t N>
 void ParticleSystemData<N>::Set(const ParticleSystemData& other)
 {
     *this = other;
+    Resize(NumberOfParticles());
 }
 
 template <size_t N>
diff --git a/Tests/PythonTests/test_mpm_system_data.py b/Tests/PythonTests/test_mpm_system_data.py
new file mode 100644
index 0000000000..f6c391fff7
--- /dev/null
+++ b/Tests/PythonTests/test_mpm_system_data.py
@@ -0,0 +1,65 @@
+import numpy as np
+import pytest
+
+import pyCubbyFlow
+
+
+@pytest.mark.parametrize(
+    "data_type,resolution,spacing,origin,grid_velocity",
+    [
+        (
+            pyCubbyFlow.MPMSystemData2,
+            pyCubbyFlow.Vector2UZ(3, 3),
+            pyCubbyFlow.Vector2D(1.0, 1.0),
+            pyCubbyFlow.Vector2D(),
+            (1.0, 1.0),
+        ),
+        (
+            pyCubbyFlow.MPMSystemData3,
+            pyCubbyFlow.Vector3UZ(3, 3, 3),
+            pyCubbyFlow.Vector3D(1.0, 1.0, 1.0),
+            pyCubbyFlow.Vector3D(),
+            (1.0, 1.0, 1.0),
+        ),
+    ],
+)
+def test_mpm_system_data_api(
+    data_type, resolution, spacing, origin, grid_velocity
+):
+    data = data_type(resolution, spacing, origin, 1)
+
+    masses = np.asarray(data.particleMasses)
+    volumes = np.asarray(data.initialVolumes)
+    assert masses.shape == (1,)
+    assert volumes.shape == (1,)
+    masses[0] = 2.0
+    volumes[0] = 3.0
+    assert data.particleMasses[0] == 2.0
+    assert data.initialVolumes[0] == 3.0
+
+    states = data.deformationStates
+    assert len(states) == 1
+    states[0].elastic[0, 0] = 2.0
+    states[0].plastic[0, 0] = 3.0
+    assert states[0].elastic[0, 0] == 2.0
+    assert states[0].plastic[0, 0] == 3.0
+
+    assert data.gridMass.resolution == resolution
+    assert data.gridVelocities.resolution == resolution
+    assert data.gridVelocitiesBeforeUpdate.resolution == resolution
+
+    data.flipBlendingFactor = 0.5
+    assert data.flipBlendingFactor == 0.5
+    data.TransferFromParticlesToGrid()
+    data.gridVelocities.Fill(grid_velocity)
+    data.TransferFromGridToParticles()
+    resized_grid = resolution + resolution, spacing + spacing, origin + spacing
+    data.ResizeGrid(*resized_grid)
+    for grid in (
+        data.gridMass,
+        data.gridVelocities,
+        data.gridVelocitiesBeforeUpdate,
+    ):
+        assert grid.resolution == resized_grid[0]
+        assert grid.gridSpacing == resized_grid[1]
+        assert grid.gridOrigin == resized_grid[2]
diff --git a/Tests/UnitTests/MPMSystemDataTests.cpp b/Tests/UnitTests/MPMSystemDataTests.cpp
new file mode 100644
index 0000000000..490e9800bb
--- /dev/null
+++ b/Tests/UnitTests/MPMSystemDataTests.cpp
@@ -0,0 +1,426 @@
+// 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/MPMSystemData.hpp>
+
+#include <limits>
+
+using namespace CubbyFlow;
+
+namespace
+{
+template <size_t N>
+void DirtyMPMState(MPMSystemData<N>& data)
+{
+    data.ParticleMasses()[0] = 7.0;
+    data.InitialVolumes()[0] = 8.0;
+    data.DeformationStates()[0].elastic(0, 0) = 9.0;
+    data.GridMass().Fill(10.0, ExecutionPolicy::Serial);
+    data.GridVelocities().Fill(Vector<double, N>::MakeConstant(11.0),
+                               ExecutionPolicy::Serial);
+    data.GridVelocitiesBeforeUpdate().Fill(
+        Vector<double, N>::MakeConstant(12.0), ExecutionPolicy::Serial);
+}
+
+template <size_t N>
+void ExpectMPMStateReset(const MPMSystemData<N>& data)
+{
+    const Vector<size_t, N> zeroIndex{};
+    const Vector<double, N> zero{};
+    const auto identity = Matrix<double, N, N>::MakeIdentity();
+
+    EXPECT_DOUBLE_EQ(data.ParticleMasses()[0], data.Mass());
+    EXPECT_DOUBLE_EQ(data.InitialVolumes()[0], 0.0);
+    EXPECT_TRUE(data.DeformationStates()[0].elastic.IsSimilar(identity));
+    EXPECT_DOUBLE_EQ(data.GridMass()(zeroIndex), 0.0);
+    EXPECT_TRUE(data.GridVelocities()(zeroIndex).IsSimilar(zero));
+    EXPECT_TRUE(data.GridVelocitiesBeforeUpdate()(zeroIndex).IsSimilar(zero));
+}
+
+template <size_t N>
+void ExpectParticleStateResizes()
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           1 };
+
+    data.ParticleMasses()[0] = 1.25;
+    data.InitialVolumes()[0] = 0.5;
+    data.SetMass(2.5);
+    data.AddParticle(Vector<double, N>::MakeConstant(0.25));
+    ASSERT_EQ(data.NumberOfParticles(), 2u);
+    EXPECT_DOUBLE_EQ(data.ParticleMasses()[0], 1.25);
+    EXPECT_DOUBLE_EQ(data.ParticleMasses()[1], 2.5);
+    EXPECT_DOUBLE_EQ(data.InitialVolumes()[0], 0.5);
+    EXPECT_DOUBLE_EQ(data.InitialVolumes()[1], 0.0);
+
+    const auto identity = Matrix<double, N, N>::MakeIdentity();
+    EXPECT_TRUE(data.DeformationStates()[0].elastic.IsSimilar(identity));
+    EXPECT_TRUE(data.DeformationStates()[1].elastic.IsSimilar(identity));
+    EXPECT_TRUE(data.DeformationStates()[0].plastic.IsSimilar(identity));
+    EXPECT_TRUE(data.DeformationStates()[1].plastic.IsSimilar(identity));
+}
+
+template <size_t N>
+void ExpectBaseSetResizesMPMState()
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(2),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           1 };
+    ParticleSystemData<N> source{ 3 };
+    source.SetMass(4.0);
+    DirtyMPMState(data);
+
+    data.Set(source);
+
+    EXPECT_EQ(data.NumberOfParticles(), 3u);
+    EXPECT_EQ(data.ParticleMasses().Length(), 3u);
+    EXPECT_EQ(data.InitialVolumes().Length(), 3u);
+    EXPECT_EQ(data.DeformationStates().Length(), 3u);
+    ExpectMPMStateReset(data);
+}
+
+template <size_t N>
+void ExpectBaseDeserializeResizesMPMState()
+{
+    ParticleSystemData<N> source{ 3 };
+    source.SetMass(4.0);
+    std::vector<uint8_t> buffer;
+    source.Serialize(&buffer);
+
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(2),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           1 };
+    DirtyMPMState(data);
+    data.Deserialize(buffer);
+
+    EXPECT_EQ(data.NumberOfParticles(), 3u);
+    EXPECT_EQ(data.ParticleMasses().Length(), 3u);
+    EXPECT_EQ(data.InitialVolumes().Length(), 3u);
+    EXPECT_EQ(data.DeformationStates().Length(), 3u);
+    ExpectMPMStateReset(data);
+}
+
+template <size_t N>
+void ExpectGridStateResizes()
+{
+    MPMSystemData<N> data;
+    const auto resolution = Vector<size_t, N>::MakeConstant(3);
+    const auto spacing = Vector<double, N>::MakeConstant(0.5);
+    const auto origin = Vector<double, N>::MakeConstant(-1.0);
+
+    data.ResizeGrid(resolution, spacing, origin);
+
+    for (size_t axis = 0; axis < N; ++axis)
+    {
+        EXPECT_EQ(data.GridMass().Resolution()[axis], resolution[axis]);
+        EXPECT_EQ(data.GridVelocities().Resolution()[axis], resolution[axis]);
+        EXPECT_EQ(data.GridVelocitiesBeforeUpdate().Resolution()[axis],
+                  resolution[axis]);
+        EXPECT_DOUBLE_EQ(data.GridMass().GridSpacing()[axis], spacing[axis]);
+        EXPECT_DOUBLE_EQ(data.GridMass().Origin()[axis], origin[axis]);
+    }
+}
+
+template <size_t N>
+void ExpectStencilPartitionAndGradient()
+{
+    const Vector<double, N> spacing = Vector<double, N>::MakeConstant(0.5);
+    const Vector<double, N> position = Vector<double, N>::MakeConstant(1.125);
+    const auto stencil = CubicBSplineKernel<N>::GetStencil(position, spacing,
+                                                           Vector<double, N>{});
+
+    double weightSum = 0.0;
+    Vector<double, N> gradientSum;
+
+    for (const auto& entry : stencil)
+    {
+        weightSum += entry.weight;
+        gradientSum += entry.gradient;
+    }
+
+    EXPECT_NEAR(weightSum, 1.0, 1e-12);
+    EXPECT_NEAR(gradientSum.Length(), 0.0, 1e-12);
+
+    constexpr double epsilon = 1e-6;
+    auto shifted = position;
+
+    shifted[0] += epsilon;
+
+    const auto shiftedStencil = CubicBSplineKernel<N>::GetStencil(
+        shifted, spacing, Vector<double, N>{});
+
+    for (size_t i = 0; i < stencil.size(); ++i)
+    {
+        EXPECT_NEAR((shiftedStencil[i].weight - stencil[i].weight) / epsilon,
+                    stencil[i].gradient[0], 1e-5);
+    }
+}
+
+template <size_t N>
+void ExpectStencilRejectsUnrepresentableCoordinates()
+{
+    const auto zero = Vector<double, N>{};
+    const auto unit = Vector<double, N>::MakeConstant(1.0);
+    const auto largest =
+        Vector<double, N>::MakeConstant(std::numeric_limits<double>::max());
+    const auto lowest =
+        Vector<double, N>::MakeConstant(std::numeric_limits<double>::lowest());
+
+    EXPECT_THROW((void)CubicBSplineKernel<N>::GetStencil(unit, zero, {}),
+                 std::invalid_argument);
+    EXPECT_THROW((void)CubicBSplineKernel<N>::GetStencil(largest, unit, {}),
+                 std::invalid_argument);
+    EXPECT_THROW((void)CubicBSplineKernel<N>::GetStencil(largest, unit, lowest),
+                 std::invalid_argument);
+}
+
+template <size_t N>
+void ExpectParticleToGridConservation(const Vector<double, N>& firstPosition)
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           2 };
+    data.Positions()[0] = firstPosition;
+    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 gridMass = 0.0;
+    Vector<double, N> gridMomentum;
+
+    data.GridMass().ForEachDataPointIndex(
+        [&data, &gridMass, &gridMomentum](const Vector<size_t, N>& index) {
+            const double mass = data.GridMass()(index);
+            gridMass += mass;
+            gridMomentum += mass * data.GridVelocities()(index);
+            EXPECT_TRUE(data.GridVelocities()(index).IsSimilar(
+                data.GridVelocitiesBeforeUpdate()(index), 1e-12));
+        });
+
+    EXPECT_NEAR(gridMass, 5.0, 1e-11);
+    EXPECT_TRUE(
+        gridMomentum.IsSimilar(Vector<double, N>::MakeConstant(1.0), 1e-11));
+}
+
+template <size_t N>
+void ExpectGridToParticleBlend(double factor, double expected,
+                               const Vector<double, N>& position,
+                               bool useDefault = false)
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(4),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           1 };
+    data.Positions()[0] = position;
+    data.Velocities()[0] = Vector<double, N>::MakeConstant(10.0);
+
+    data.GridVelocitiesBeforeUpdate().Fill(Vector<double, N>::MakeConstant(1.0),
+                                           ExecutionPolicy::Serial);
+    data.GridVelocities().Fill(Vector<double, N>::MakeConstant(3.0),
+                               ExecutionPolicy::Serial);
+
+    if (!useDefault)
+    {
+        data.SetFLIPBlendingFactor(factor);
+    }
+
+    data.TransferFromGridToParticles();
+
+    EXPECT_TRUE(data.Velocities()[0].IsSimilar(
+        Vector<double, N>::MakeConstant(expected), 1e-12));
+}
+
+template <size_t N>
+void ExpectRejectsDivergentGridState()
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(4) };
+    data.GridMass().Clear();
+    EXPECT_THROW(data.TransferFromParticlesToGrid(), std::invalid_argument);
+
+    data.ResizeGrid(Vector<size_t, N>::MakeConstant(4),
+                    Vector<double, N>::MakeConstant(1.0), {});
+    data.GridVelocitiesBeforeUpdate().Resize(
+        Vector<size_t, N>::MakeConstant(2),
+        Vector<double, N>::MakeConstant(1.0), {});
+    EXPECT_THROW(data.TransferFromGridToParticles(), std::invalid_argument);
+}
+
+template <size_t N>
+void ExpectInvalidGridLeavesParticleVelocitiesUnchanged()
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(8),
+                           Vector<double, N>::MakeConstant(1.0),
+                           {},
+                           2 };
+    data.Positions()[0] = Vector<double, N>::MakeConstant(1.0);
+    data.Positions()[1] = Vector<double, N>::MakeConstant(7.0);
+    data.Velocities()[0] = Vector<double, N>::MakeConstant(10.0);
+    data.Velocities()[1] = Vector<double, N>::MakeConstant(20.0);
+    data.GridVelocitiesBeforeUpdate().Fill(Vector<double, N>::MakeConstant(1.0),
+                                           ExecutionPolicy::Serial);
+    data.GridVelocities().Fill(Vector<double, N>::MakeConstant(2.0),
+                               ExecutionPolicy::Serial);
+
+    auto invalidVelocity = Vector<double, N>{};
+    invalidVelocity[0] = std::numeric_limits<double>::infinity();
+    data.GridVelocities()(Vector<size_t, N>::MakeConstant(7)) = invalidVelocity;
+
+    EXPECT_THROW(data.TransferFromGridToParticles(), std::invalid_argument);
+    EXPECT_TRUE(data.Velocities()[0].IsSimilar(
+        Vector<double, N>::MakeConstant(10.0), 1e-12));
+    EXPECT_TRUE(data.Velocities()[1].IsSimilar(
+        Vector<double, N>::MakeConstant(20.0), 1e-12));
+}
+
+template <size_t N>
+void ExpectRejectsInvalidInput()
+{
+    MPMSystemData<N> data{ Vector<size_t, N>::MakeConstant(2) };
+    const Vector<size_t, N> zeroResolution{};
+    const auto unitSpacing = Vector<double, N>::MakeConstant(1.0);
+    auto negativeSpacing = unitSpacing;
+    negativeSpacing[0] = -1.0;
+    const auto overflowingResolution =
+        Vector<size_t, N>::MakeConstant(std::numeric_limits<size_t>::max());
+
+    EXPECT_THROW(data.ResizeGrid(zeroResolution, unitSpacing, {}),
+                 std::invalid_argument);
+    EXPECT_THROW(data.ResizeGrid(Vector<size_t, N>::MakeConstant(2),
+                                 negativeSpacing, {}),
+                 std::invalid_argument);
+    EXPECT_THROW(data.ResizeGrid(overflowingResolution, unitSpacing, {}),
+                 std::invalid_argument);
+    EXPECT_THROW(data.SetFLIPBlendingFactor(-0.1), std::invalid_argument);
+    EXPECT_THROW(data.SetFLIPBlendingFactor(1.1), std::invalid_argument);
+    EXPECT_THROW(
+        data.SetFLIPBlendingFactor(std::numeric_limits<double>::quiet_NaN()),
+        std::invalid_argument);
+
+    data.Resize(1);
+    data.ParticleMasses()[0] = 0.0;
+    EXPECT_THROW(data.TransferFromParticlesToGrid(), std::invalid_argument);
+    data.ParticleMasses()[0] = 1.0;
+    data.Positions()[0][0] = std::numeric_limits<double>::infinity();
+    EXPECT_THROW(data.TransferFromParticlesToGrid(), std::invalid_argument);
+    EXPECT_THROW(data.TransferFromGridToParticles(), std::invalid_argument);
+
+    data.Positions()[0] = Vector<double, N>{};
+    data.Velocities()[0] = Vector<double, N>{};
+    auto invalidGridVelocity = Vector<double, N>{};
+    invalidGridVelocity[0] = std::numeric_limits<double>::infinity();
+    data.GridVelocities().Fill(invalidGridVelocity, ExecutionPolicy::Serial);
+    EXPECT_THROW(data.TransferFromGridToParticles(), std::invalid_argument);
+}
+}  // namespace
+
+TEST(CubicBSplineKernel, Values)
+{
+    EXPECT_NEAR(CubicBSplineKernel<2>::Weight(0.0), 2.0 / 3.0, 1e-12);
+    EXPECT_NEAR(CubicBSplineKernel<2>::Weight(1.0), 1.0 / 6.0, 1e-12);
+    EXPECT_DOUBLE_EQ(CubicBSplineKernel<2>::Weight(2.0), 0.0);
+    EXPECT_NEAR(CubicBSplineKernel<2>::Gradient(0.5), -0.625, 1e-12);
+    EXPECT_NEAR(CubicBSplineKernel<2>::Gradient(-0.5), 0.625, 1e-12);
+}
+
+TEST(CubicBSplineKernel, Stencil)
+{
+    ExpectStencilPartitionAndGradient<2>();
+    ExpectStencilPartitionAndGradient<3>();
+}
+
+TEST(CubicBSplineKernel, RejectsUnrepresentableCoordinates)
+{
+    ExpectStencilRejectsUnrepresentableCoordinates<2>();
+    ExpectStencilRejectsUnrepresentableCoordinates<3>();
+}
+
+TEST(MPMSystemData, ParticleStateResizes)
+{
+    ExpectParticleStateResizes<2>();
+    ExpectParticleStateResizes<3>();
+}
+
+TEST(MPMSystemData, BaseSetResizesMPMState)
+{
+    ExpectBaseSetResizesMPMState<2>();
+    ExpectBaseSetResizesMPMState<3>();
+}
+
+TEST(MPMSystemData, BaseDeserializeResizesMPMState)
+{
+    ExpectBaseDeserializeResizesMPMState<2>();
+    ExpectBaseDeserializeResizesMPMState<3>();
+}
+
+TEST(MPMSystemData, GridStateResizes)
+{
+    ExpectGridStateResizes<2>();
+    ExpectGridStateResizes<3>();
+}
+
+TEST(MPMSystemData, ParticleToGridConservesMassAndMomentum)
+{
+    ExpectParticleToGridConservation<2>(Vector2D::MakeConstant(1.25));
+    ExpectParticleToGridConservation<3>(Vector3D::MakeConstant(1.25));
+    ExpectParticleToGridConservation<2>(Vector2D{});
+    ExpectParticleToGridConservation<3>(Vector3D{});
+}
+
+TEST(MPMSystemData, GridToParticleBlendsPICAndFLIP)
+{
+    ExpectGridToParticleBlend<2>(0.0, 3.0, Vector2D::MakeConstant(1.25));
+    ExpectGridToParticleBlend<3>(0.0, 3.0, Vector3D::MakeConstant(1.25));
+    ExpectGridToParticleBlend<2>(1.0, 12.0, Vector2D{});
+    ExpectGridToParticleBlend<3>(1.0, 12.0, Vector3D{});
+    ExpectGridToParticleBlend<2>(0.95, 11.55, Vector2D::MakeConstant(1.25),
+                                 true);
+    ExpectGridToParticleBlend<3>(0.95, 11.55, Vector3D::MakeConstant(1.25),
+                                 true);
+}
+
+TEST(MPMSystemData, EmptyTransfersAreSafe)
+{
+    MPMSystemData2 data2;
+    MPMSystemData3 data3;
+
+    EXPECT_NO_THROW(data2.TransferFromParticlesToGrid());
+    EXPECT_NO_THROW(data2.TransferFromGridToParticles());
+    EXPECT_NO_THROW(data3.TransferFromParticlesToGrid());
+    EXPECT_NO_THROW(data3.TransferFromGridToParticles());
+}
+
+TEST(MPMSystemData, RejectsDivergentGridState)
+{
+    ExpectRejectsDivergentGridState<2>();
+    ExpectRejectsDivergentGridState<3>();
+}
+
+TEST(MPMSystemData, InvalidGridLeavesParticleVelocitiesUnchanged)
+{
+    ExpectInvalidGridLeavesParticleVelocitiesUnchanged<2>();
+    ExpectInvalidGridLeavesParticleVelocitiesUnchanged<3>();
+}
+
+TEST(MPMSystemData, RejectsInvalidInput)
+{
+    ExpectRejectsInvalidInput<2>();
+    ExpectRejectsInvalidInput<3>();
+}