Implement explicit 2-D and 3-D snow MPM solvers

Includes/API/Python/Solver/Particle/MPM/SnowMPMSolver.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_SNOW_MPM_SOLVER_HPP
+#define CUBBYFLOW_PYTHON_SNOW_MPM_SOLVER_HPP
+
+#include <pybind11/pybind11.h>
+
+void AddSnowMPMSolver2(pybind11::module& m);
+void AddSnowMPMSolver3(pybind11::module& m);
+
+#endif

Includes/Core/Particle/MPM/SnowConstitutiveModel-Impl.hpp

@@ -113,13 +113,7 @@ SnowConstitutiveModel<N>::ComputeKirchhoffStress(const State& state) const
 
     const MatrixType rotation = u * v.Transposed();
     const double elasticDeterminant = state.elastic.Determinant();
-    const double hardening =
-        std::exp(m_hardeningCoefficient * (1.0 - state.plastic.Determinant()));
-
-    if (!std::isfinite(hardening))
-    {
-        throw std::invalid_argument{ "Non-finite snow hardening." };
-    }
+    const double hardening = ComputeHardening(state);
 
     const double mu = m_mu0 * hardening;
     const double lambda = m_lambda0 * hardening;
@@ -136,6 +130,80 @@ SnowConstitutiveModel<N>::ComputeKirchhoffStress(const State& state) const
     return stress;
 }
 
+template <size_t N>
+double SnowConstitutiveModel<N>::ComputeWaveSpeed(const State& state,
+                                                  double referenceDensity) const
+{
+    ValidateDeformation(state.elastic);
+    ValidateDeformation(state.plastic);
+
+    if (!std::isfinite(referenceDensity) || referenceDensity <= 0.0)
+    {
+        throw std::invalid_argument{ "Invalid snow reference density." };
+    }
+
+    MatrixType u;
+    Vector<double, N> singularValues;
+    MatrixType v;
+    SVD(state.elastic, u, singularValues, v);
+
+    const double elasticDeterminant = state.elastic.Determinant();
+    const double totalDeterminant =
+        elasticDeterminant * state.plastic.Determinant();
+    const double currentDensity = referenceDensity / totalDeterminant;
+    const double hardening = ComputeHardening(state);
+    const double mu = m_mu0 * hardening;
+    const double lambda = m_lambda0 * hardening;
+    double maxCandidate = 0.0;
+
+    for (size_t a = 0; a < N; ++a)
+    {
+        const double sigmaA = singularValues[a];
+        maxCandidate = std::max(
+            maxCandidate, 2.0 * mu * sigmaA * sigmaA +
+                              lambda * elasticDeterminant * elasticDeterminant);
+
+        for (size_t b = 0; b < N; ++b)
+        {
+            if (a == b)
+            {
+                continue;
+            }
+
+            const double sigmaB = singularValues[b];
+            maxCandidate = std::max(maxCandidate, 2.0 * mu * sigmaB * sigmaB *
+                                                      (sigmaA + sigmaB - 1.0) /
+                                                      (sigmaA + sigmaB));
+        }
+    }
+
+    const double kappa = maxCandidate / totalDeterminant;
+    const double waveSpeed = std::sqrt(kappa / currentDensity);
+
+    if (!std::isfinite(totalDeterminant) || totalDeterminant <= 0.0 ||
+        !std::isfinite(currentDensity) || currentDensity <= 0.0 ||
+        !std::isfinite(kappa) || kappa <= 0.0 || !std::isfinite(waveSpeed))
+    {
+        throw std::invalid_argument{ "Invalid snow wave speed." };
+    }
+
+    return waveSpeed;
+}
+
+template <size_t N>
+double SnowConstitutiveModel<N>::ComputeHardening(const State& state) const
+{
+    const double hardening =
+        std::exp(m_hardeningCoefficient * (1.0 - state.plastic.Determinant()));
+
+    if (!std::isfinite(hardening))
+    {
+        throw std::invalid_argument{ "Non-finite snow hardening." };
+    }
+
+    return hardening;
+}
+
 template <size_t N>
 bool SnowConstitutiveModel<N>::IsFinite(const MatrixType& matrix)
 {

Includes/Core/Particle/MPM/SnowConstitutiveModel.hpp

@@ -90,7 +90,22 @@ class SnowConstitutiveModel final
     //! \return Fixed-corotated Kirchhoff stress.
     [[nodiscard]] MatrixType ComputeKirchhoffStress(const State& state) const;
 
+    //!
+    //! \brief Estimates the fastest elastic wave speed for the state.
+    //!
+    //! Evaluates the principal-basis acoustic-tensor candidates for the
+    //! fixed-corotated snow energy.
+    //!
+    //! \param[in] state Current deformation state.
+    //! \param[in] referenceDensity Positive material density at rest.
+    //!
+    //! \return Estimated elastic wave speed.
+    [[nodiscard]] double ComputeWaveSpeed(const State& state,
+                                          double referenceDensity) const;
+
  private:
+    [[nodiscard]] double ComputeHardening(const State& state) const;
+
     [[nodiscard]] static bool IsFinite(const MatrixType& matrix);
 
     static void ValidateDeformation(const MatrixType& matrix);