diff --git a/Includes/Core/Solver/Particle/MPM/MPMFluidSolver-Impl.hpp b/Includes/Core/Solver/Particle/MPM/MPMFluidSolver-Impl.hpp index d90cbe854e..74b482f979 100644 --- a/Includes/Core/Solver/Particle/MPM/MPMFluidSolver-Impl.hpp +++ b/Includes/Core/Solver/Particle/MPM/MPMFluidSolver-Impl.hpp @@ -12,6 +12,7 @@ #define CUBBYFLOW_MPM_FLUID_SOLVER_IMPL_HPP #include +#include #include #include #include @@ -76,6 +77,18 @@ void MPMFluidSolver::SetTimeStepLimitScale(double newScale) m_timeStepLimitScale = newScale; } +template +int MPMFluidSolver::GetClosedDomainBoundaryFlag() const +{ + return m_closedDomainBoundaryFlag; +} + +template +void MPMFluidSolver::SetClosedDomainBoundaryFlag(int flag) +{ + m_closedDomainBoundaryFlag = flag; +} + template MPMFluidSolver::Builder MPMFluidSolver::GetBuilder() { @@ -154,11 +167,21 @@ template void MPMFluidSolver::OnBeginAdvanceTimeStep(double timeStepInSeconds) { m_mpmSystemData->TransferFromParticlesToGrid(); + InitializeReferenceVolumes(); UpdateGridVelocities(timeStepInSeconds); + ConstrainGridVelocities(); + m_mpmSystemData->TransferFromGridToParticles(timeStepInSeconds); } +template +void MPMFluidSolver::OnEndAdvanceTimeStep(double timeStepInSeconds) +{ + Base::OnEndAdvanceTimeStep(timeStepInSeconds); + ConstrainParticlesToDomain(); +} + template MPMFluidSolver::SizeType MPMFluidSolver::ClampIndex( const Vector& index, const SizeType& dataSize) @@ -252,6 +275,86 @@ void MPMFluidSolver::UpdateGridVelocities(double timeStepInSeconds) } } +template +void MPMFluidSolver::ConstrainGridVelocities() +{ + static constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, + DIRECTION_BACK }; + static constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP, + DIRECTION_FRONT }; + + const auto& gridMass = m_mpmSystemData->GridMass(); + auto& gridVelocities = m_mpmSystemData->GridVelocities(); + const auto dataSize = gridVelocities.DataSize(); + + gridVelocities.ParallelForEachDataPointIndex( + [this, &gridMass, &gridVelocities, dataSize](const SizeType& index) { + if (gridMass(index) <= 0.0) + { + return; + } + + VectorType velocity = gridVelocities(index); + + if (const auto& collider = this->GetCollider(); collider != nullptr) + { + VectorType position = gridVelocities.DataPosition()(index); + collider->ResolveCollision(0.0, 0.0, &position, &velocity); + } + + for (size_t axis = 0; axis < N; ++axis) + { + const bool exceedsLower = + (m_closedDomainBoundaryFlag & lowerFlags[axis]) != 0 && + index[axis] == 0 && velocity[axis] < 0.0; + const bool exceedsUpper = + (m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 && + index[axis] == dataSize[axis] - 1 && velocity[axis] > 0.0; + + if (exceedsLower || exceedsUpper) + { + velocity[axis] = 0.0; + } + } + + gridVelocities(index) = velocity; + }); +} + +template +void MPMFluidSolver::ConstrainParticlesToDomain() +{ + static constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, + DIRECTION_BACK }; + static constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP, + DIRECTION_FRONT }; + + const auto domain = m_mpmSystemData->GridMass().GetBoundingBox(); + auto positions = m_mpmSystemData->Positions(); + auto velocities = m_mpmSystemData->Velocities(); + + ParallelFor( + ZERO_SIZE, positions.Length(), + [&domain, &positions, &velocities, this](size_t i) { + for (size_t axis = 0; axis < N; ++axis) + { + if ((m_closedDomainBoundaryFlag & lowerFlags[axis]) != 0 && + positions[i][axis] <= domain.lowerCorner[axis]) + { + positions[i][axis] = domain.lowerCorner[axis]; + velocities[i][axis] = std::max(velocities[i][axis], 0.0); + } + + if ((m_closedDomainBoundaryFlag & upperFlags[axis]) != 0 && + positions[i][axis] >= domain.upperCorner[axis]) + { + positions[i][axis] = domain.upperCorner[axis]; + velocities[i][axis] = std::min(velocities[i][axis], 0.0); + } + } + }); +} + template MPMFluidSolver::Builder& MPMFluidSolver::Builder::WithResolution( const SizeType& resolution) diff --git a/Includes/Core/Solver/Particle/MPM/MPMFluidSolver.hpp b/Includes/Core/Solver/Particle/MPM/MPMFluidSolver.hpp index af34904d93..114b92d5ed 100644 --- a/Includes/Core/Solver/Particle/MPM/MPMFluidSolver.hpp +++ b/Includes/Core/Solver/Particle/MPM/MPMFluidSolver.hpp @@ -28,6 +28,13 @@ namespace CubbyFlow //! Uses the Tait equation of state with existing MPM particle-grid transfers //! and particle-solver collision lifecycle. //! +//! \note Choose the artificial speed of sound high enough to keep the maximum +//! relative density error near or below 3%; ten times the expected peak flow +//! speed is a common starting point. Adaptive stepping limits each step to +//! `timeStepLimitScale * minGridSpacing / (speedOfSound + maxParticleSpeed)`. +//! This explicit limit assumes weak compression; reduce the scale or increase +//! the sound speed when density error or pressure oscillation is excessive. +//! template class MPMFluidSolver : public std::conditional_t @@ -67,6 +74,12 @@ class MPMFluidSolver : public std::conditional_t& index, const SizeType& dataSize); @@ -92,9 +108,14 @@ class MPMFluidSolver : public std::conditional_t> m_mpmSystemData; MPMFluidConstitutiveModel m_constitutiveModel; double m_timeStepLimitScale = 0.9; + int m_closedDomainBoundaryFlag = DIRECTION_ALL; }; //! Front-end to create MPMFluidSolver objects step by step. diff --git a/Sources/API/Python/Solver/Particle/MPM/MPMFluidSolver.cpp b/Sources/API/Python/Solver/Particle/MPM/MPMFluidSolver.cpp index 91c888bea8..3ffccfebe6 100644 --- a/Sources/API/Python/Solver/Particle/MPM/MPMFluidSolver.cpp +++ b/Sources/API/Python/Solver/Particle/MPM/MPMFluidSolver.cpp @@ -38,6 +38,9 @@ void AddMPMFluidSolver(pybind11::module& m, const char* name) .def_property_readonly("mpmSystemData", &Solver::GetMPMSystemData) .def_property("timeStepLimitScale", &Solver::GetTimeStepLimitScale, &Solver::SetTimeStepLimitScale) + .def_property("closedDomainBoundaryFlag", + &Solver::GetClosedDomainBoundaryFlag, + &Solver::SetClosedDomainBoundaryFlag) .def_property_readonly( "targetDensity", [](const Solver& solver) { diff --git a/Tests/ManualTests/MPMFluidSolverTests.cpp b/Tests/ManualTests/MPMFluidSolverTests.cpp new file mode 100644 index 0000000000..3b7fa27060 --- /dev/null +++ b/Tests/ManualTests/MPMFluidSolverTests.cpp @@ -0,0 +1,79 @@ +// 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 + +#include +#include +#include + +using namespace CubbyFlow; + +TEST(MPMFluidSolverManual, DamBreakReportsDensityErrorAndMomentum) +{ + constexpr double targetDensity = 1000.0; + constexpr double particleSpacing = 0.05; + constexpr double particleMass = + targetDensity * particleSpacing * particleSpacing; + MPMFluidSolver2 solver{ { 10, 10 }, { 0.1, 0.1 }, {}, 0.02, + particleMass, targetDensity, 30.0 }; + solver.SetDragCoefficient(0.0); + solver.SetTimeStepLimitScale(0.5); + + Array1 positions; + for (size_t j = 0; j < 6; ++j) + { + for (size_t i = 0; i < 3; ++i) + { + positions.Append( + { 0.05 + particleSpacing * static_cast(i), + 0.05 + particleSpacing * static_cast(j) }); + } + } + + auto data = solver.GetMPMSystemData(); + data->AddParticles(positions); + + for (Frame frame{ 0, 1.0 / 240.0 }; frame.index < 40; ++frame) + { + solver.Update(frame); + } + + const auto masses = data->ParticleMasses(); + const auto volumes = data->InitialVolumes(); + const auto volumeRatios = data->VolumeRatios(); + const auto velocities = data->Velocities(); + double maxDensityError = 0.0; + Vector2D totalMomentum; + + for (size_t i = 0; i < data->NumberOfParticles(); ++i) + { + const double density = solver.GetConstitutiveModel().ComputeDensity( + masses[i], volumes[i] * volumeRatios[i]); + maxDensityError = + std::max(maxDensityError, std::abs(density / targetDensity - 1.0)); + totalMomentum += masses[i] * velocities[i]; + + EXPECT_TRUE(std::isfinite(density)); + EXPECT_TRUE(std::isfinite(velocities[i].x)); + EXPECT_TRUE(std::isfinite(velocities[i].y)); + } + + std::cout << "MPM dam-break max density error: " << 100.0 * maxDensityError + << "%\n" + << "MPM dam-break total momentum: (" << totalMomentum.x << ", " + << totalMomentum.y << ")\n"; + + EXPECT_TRUE(std::isfinite(maxDensityError)); + EXPECT_TRUE(std::isfinite(totalMomentum.x)); + EXPECT_TRUE(std::isfinite(totalMomentum.y)); +} diff --git a/Tests/PythonTests/test_mpm_fluid_solver.py b/Tests/PythonTests/test_mpm_fluid_solver.py index 40d02f6397..efbd26afb6 100644 --- a/Tests/PythonTests/test_mpm_fluid_solver.py +++ b/Tests/PythonTests/test_mpm_fluid_solver.py @@ -47,9 +47,12 @@ def test_mpm_fluid_solver_api( assert solver.eosExponent == pytest.approx(5.0) assert solver.negativePressureScale == pytest.approx(0.25) assert solver.timeStepLimitScale == pytest.approx(0.9) + assert solver.closedDomainBoundaryFlag == pyCubbyFlow.DIRECTION_ALL solver.timeStepLimitScale = 0.5 assert solver.timeStepLimitScale == pytest.approx(0.5) + solver.closedDomainBoundaryFlag = pyCubbyFlow.DIRECTION_LEFT + assert solver.closedDomainBoundaryFlag == pyCubbyFlow.DIRECTION_LEFT with pytest.raises(AttributeError): solver.targetDensity = 1000.0 diff --git a/Tests/UnitTests/MPMFluidSolverTests.cpp b/Tests/UnitTests/MPMFluidSolverTests.cpp index 0d80b5a34c..f37a158f98 100644 --- a/Tests/UnitTests/MPMFluidSolverTests.cpp +++ b/Tests/UnitTests/MPMFluidSolverTests.cpp @@ -12,9 +12,13 @@ #include #include +#include +#include #include #include +#include +#include #include #include #include @@ -67,6 +71,11 @@ template void ExpectParameters() { MPMFluidSolver solver; + EXPECT_EQ(solver.GetClosedDomainBoundaryFlag(), DIRECTION_ALL); + + solver.SetClosedDomainBoundaryFlag(DIRECTION_LEFT | DIRECTION_UP); + EXPECT_EQ(solver.GetClosedDomainBoundaryFlag(), + DIRECTION_LEFT | DIRECTION_UP); EXPECT_DOUBLE_EQ(solver.GetTimeStepLimitScale(), 0.9); EXPECT_DOUBLE_EQ(solver.GetConstitutiveModel().GetTargetDensity(), WATER_DENSITY); @@ -333,6 +342,390 @@ void ExpectCompressedParticlesMoveOutwardAndStayFinite() } } } + +template +void ExpectClosedDomainWalls() +{ + constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, + DIRECTION_BACK }; + constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP, + DIRECTION_FRONT }; + const auto resolution = VectorUZ::MakeConstant(4); + const auto spacing = VectorD::MakeConstant(1.0); + + for (size_t axis = 0; axis < N; ++axis) + { + for (bool isUpper : { false, true }) + { + TestableMPMFluidSolver solver{ resolution, spacing }; + solver.SetClosedDomainBoundaryFlag(isUpper ? upperFlags[axis] + : lowerFlags[axis]); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD position = VectorD::MakeConstant(2.0); + VectorD velocity; + + position[axis] = isUpper ? 3.75 : 0.25; + velocity[axis] = isUpper ? 1.0 : -1.0; + + auto data = solver.GetMPMSystemData(); + data->AddParticle(position, velocity); + + solver.BeginStep(1e-3); + + VectorUZ nodeIndex = VectorUZ::MakeConstant(2); + nodeIndex[axis] = + isUpper ? data->GridMass().DataSize()[axis] - 1 : 0; + + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); + EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[axis], 0.0); + } + } + + TestableMPMFluidSolver open{ resolution, spacing }; + open.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + open.SetGravity({}); + open.SetDragCoefficient(0.0); + + VectorD position = VectorD::MakeConstant(2.0); + VectorD velocity; + + position[0] = 0.25; + velocity[0] = -1.0; + + auto data = open.GetMPMSystemData(); + data->AddParticle(position, velocity); + + open.BeginStep(1e-3); + + VectorUZ nodeIndex = VectorUZ::MakeConstant(2); + nodeIndex[0] = 0; + + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); + EXPECT_NEAR(data->GridVelocities()(nodeIndex)[0], -1.0, 1e-12); +} + +template +void ExpectMovingColliderAffectsGrid() +{ + TestableMPMFluidSolver solver{ VectorUZ::MakeConstant(4), + VectorD::MakeConstant(1.0) }; + solver.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD normal; + normal[0] = 1.0; + + auto collider = std::make_shared>( + std::make_shared>(normal, VectorD{})); + collider->linearVelocity[0] = 1.0; + + solver.SetCollider(collider); + + VectorD position = VectorD::MakeConstant(2.0); + position[0] = 0.25; + + auto data = solver.GetMPMSystemData(); + data->AddParticle(position); + + solver.BeginStep(1e-3); + + VectorUZ nodeIndex = VectorUZ::MakeConstant(2); + nodeIndex[0] = 0; + + ASSERT_GT(data->GridMass()(nodeIndex), 0.0); + EXPECT_DOUBLE_EQ(data->GridVelocities()(nodeIndex)[0], 1.0); +} + +template +void ExpectFrictionAffectsGrid() +{ + const auto runCase = [](double frictionCoefficient) { + TestableMPMFluidSolver solver{ VectorUZ::MakeConstant(4), + VectorD::MakeConstant(1.0) }; + solver.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD normal; + normal[1] = 1.0; + + VectorD point; + point[1] = 2.0; + + auto collider = std::make_shared>( + std::make_shared>(normal, point)); + collider->SetFrictionCoefficient(frictionCoefficient); + + solver.SetCollider(collider); + + VectorD position = VectorD::MakeConstant(2.25); + VectorD velocity; + + velocity[0] = 1.0; + velocity[1] = -1.0; + + auto data = solver.GetMPMSystemData(); + data->AddParticle(position, velocity); + + solver.BeginStep(1e-3); + + const VectorUZ nodeIndex = VectorUZ::MakeConstant(2); + EXPECT_GT(data->GridMass()(nodeIndex), 0.0); + + return data->GridVelocities()(nodeIndex); + }; + + const auto frictionless = runCase(0.0); + EXPECT_NEAR(frictionless[0], 1.0, 1e-12); + EXPECT_DOUBLE_EQ(frictionless[1], 0.0); + + const auto frictional = runCase(1.0); + EXPECT_DOUBLE_EQ(frictional[0], 0.0); + EXPECT_DOUBLE_EQ(frictional[1], 0.0); +} + +template +void ExpectParticleDomainProjection(size_t axis, bool isUpper) +{ + constexpr std::array lowerFlags{ DIRECTION_LEFT, DIRECTION_DOWN, + DIRECTION_BACK }; + constexpr std::array upperFlags{ DIRECTION_RIGHT, DIRECTION_UP, + DIRECTION_FRONT }; + const auto runCase = [axis, isUpper](int boundaryFlag, + double normalVelocity) { + MPMFluidSolver solver{ VectorUZ::MakeConstant(4), + VectorD::MakeConstant(1.0) }; + + UseOneFixedStep(&solver); + + solver.SetClosedDomainBoundaryFlag(boundaryFlag); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD position = VectorD::MakeConstant(2.0); + VectorD velocity; + + position[axis] = isUpper ? 4.01 : -0.01; + velocity[axis] = normalVelocity; + + auto data = solver.GetMPMSystemData(); + data->AddParticle(position, velocity); + + solver.Update(Frame{ 0, 1e-3 }); + + return std::array{ data->Positions()[0][axis], + data->Velocities()[0][axis] }; + }; + + const int boundaryFlag = isUpper ? upperFlags[axis] : lowerFlags[axis]; + const double outwardVelocity = isUpper ? 1.0 : -1.0; + const double boundary = isUpper ? 4.0 : 0.0; + + const auto closed = runCase(boundaryFlag, outwardVelocity); + EXPECT_DOUBLE_EQ(closed[0], boundary); + EXPECT_NEAR(closed[1], 0.0, 1e-12); + + const auto inward = runCase(boundaryFlag, -outwardVelocity); + EXPECT_DOUBLE_EQ(inward[0], boundary); + EXPECT_NEAR(inward[1], -outwardVelocity, 1e-12); + + const auto open = runCase(DIRECTION_NONE, outwardVelocity); + EXPECT_NEAR(open[0], isUpper ? 4.011 : -0.011, 1e-12); + EXPECT_NEAR(open[1], outwardVelocity, 1e-12); +} + +template +void ExpectParticleColliderProjection() +{ + constexpr double radius = 0.1; + MPMFluidSolver solver{ VectorUZ::MakeConstant(4), + VectorD::MakeConstant(1.0), + {}, + radius, + 1.0 }; + + UseOneFixedStep(&solver); + + solver.SetClosedDomainBoundaryFlag(DIRECTION_NONE); + solver.SetGravity({}); + solver.SetDragCoefficient(0.0); + + VectorD normal; + normal[1] = 1.0; + + VectorD point; + point[1] = 1.0; + + solver.SetCollider(std::make_shared>( + std::make_shared>(normal, point))); + + VectorD position = VectorD::MakeConstant(2.0); + VectorD velocity; + + position[1] = 0.95; + velocity[1] = -1.0; + + auto data = solver.GetMPMSystemData(); + data->AddParticle(position, velocity); + + solver.Update(Frame{ 0, 1e-3 }); + + EXPECT_GE(data->Positions()[0][1], point[1] + radius); + EXPECT_GE(data->Velocities()[0][1], 0.0); +} + +template +double ParticleMeasure(double spacing) +{ + double result = 1.0; + + for (size_t axis = 0; axis < N; ++axis) + { + result *= spacing; + } + + return result; +} + +template +void AddParticleBlock(MPMFluidSolver* solver, const VectorUZ& size, + double spacing) +{ + Array1> positions; + + ForEachIndex(size, [&positions, spacing](auto... rawIndices) { + const VectorUZ index{ rawIndices... }; + VectorD position; + + for (size_t axis = 0; axis < N; ++axis) + { + position[axis] = 0.05 + spacing * static_cast(index[axis]); + } + + positions.Append(position); + }); + + solver->GetMPMSystemData()->AddParticles(positions); +} + +template +void ExpectHydrostaticColumnRemainsBounded() +{ + constexpr double particleSpacing = 0.1; + constexpr double targetDensity = 1000.0; + const double particleMass = + targetDensity * ParticleMeasure(particleSpacing); + + MPMFluidSolver solver{ VectorUZ::MakeConstant(6), + VectorD::MakeConstant(0.1), + {}, + 0.02, + particleMass, + targetDensity, + 50.0 }; + solver.SetDragCoefficient(0.0); + solver.SetTimeStepLimitScale(0.5); + + VectorUZ blockSize = VectorUZ::MakeConstant(6); + blockSize[1] = 4; + + AddParticleBlock(&solver, blockSize, particleSpacing); + + for (int frame = 0; frame < 24; ++frame) + { + solver.Update(Frame{ frame, 1.0 / 240.0 }); + } + + const auto data = solver.GetMPMSystemData(); + const auto masses = data->ParticleMasses(); + const auto initialVolumes = data->InitialVolumes(); + const auto volumeRatios = data->VolumeRatios(); + double maxDensityError = 0.0; + + for (size_t i = 0; i < data->NumberOfParticles(); ++i) + { + const double density = solver.GetConstitutiveModel().ComputeDensity( + masses[i], initialVolumes[i] * volumeRatios[i]); + const double pressure = + solver.GetConstitutiveModel().ComputePressure(density); + + EXPECT_TRUE(std::isfinite(density)); + EXPECT_TRUE(std::isfinite(pressure)); + EXPECT_GT(density, 0.0); + + maxDensityError = + std::max(maxDensityError, std::abs(density / targetDensity - 1.0)); + } + + EXPECT_LT(maxDensityError, 0.03); +} + +template +void ExpectDamBreakConservesMassAndStaysInDomain() +{ + constexpr double particleSpacing = 0.05; + constexpr double targetDensity = 1000.0; + const double particleMass = + targetDensity * ParticleMeasure(particleSpacing); + + MPMFluidSolver solver{ VectorUZ::MakeConstant(10), + VectorD::MakeConstant(0.1), + {}, + 0.02, + particleMass, + targetDensity, + 30.0 }; + solver.SetDragCoefficient(0.0); + solver.SetTimeStepLimitScale(0.5); + + VectorUZ blockSize = VectorUZ::MakeConstant(3); + blockSize[1] = 6; + + AddParticleBlock(&solver, blockSize, particleSpacing); + + const auto data = solver.GetMPMSystemData(); + const size_t initialCount = data->NumberOfParticles(); + const double initialMass = particleMass * static_cast(initialCount); + double initialMaxX = 0.0; + + for (const auto& position : data->Positions()) + { + initialMaxX = std::max(initialMaxX, position[0]); + } + + for (int frame = 0; frame < 40; ++frame) + { + solver.Update(Frame{ frame, 1.0 / 240.0 }); + } + + const auto domain = data->GridMass().GetBoundingBox(); + double finalMass = 0.0; + double finalMaxX = 0.0; + + for (size_t i = 0; i < data->NumberOfParticles(); ++i) + { + finalMass += data->ParticleMasses()[i]; + finalMaxX = std::max(finalMaxX, data->Positions()[i][0]); + + for (size_t axis = 0; axis < N; ++axis) + { + EXPECT_GE(data->Positions()[i][axis], domain.lowerCorner[axis]); + EXPECT_LE(data->Positions()[i][axis], domain.upperCorner[axis]); + } + + EXPECT_TRUE(std::isfinite(data->VolumeRatios()[i])); + EXPECT_GT(data->VolumeRatios()[i], 0.0); + } + + EXPECT_EQ(data->NumberOfParticles(), initialCount); + EXPECT_NEAR(finalMass, initialMass, + std::numeric_limits::epsilon() * initialMass * + static_cast(initialCount)); + EXPECT_GT(finalMaxX, initialMaxX + 0.01); +} } // namespace TEST(MPMFluidSolver, ParametersAndBuilder) @@ -372,3 +765,50 @@ TEST(MPMFluidSolver, CompressedParticlesMoveOutwardAndStayFinite) ExpectCompressedParticlesMoveOutwardAndStayFinite<2>(); ExpectCompressedParticlesMoveOutwardAndStayFinite<3>(); } + +TEST(MPMFluidSolver, ClosedDomainWalls) +{ + ExpectClosedDomainWalls<2>(); + ExpectClosedDomainWalls<3>(); +} + +TEST(MPMFluidSolver, ColliderContact) +{ + ExpectMovingColliderAffectsGrid<2>(); + ExpectMovingColliderAffectsGrid<3>(); + ExpectFrictionAffectsGrid<2>(); + ExpectFrictionAffectsGrid<3>(); +} + +TEST(MPMFluidSolver, ParticleDomainProjection) +{ + for (bool isUpper : { false, true }) + { + for (size_t axis = 0; axis < 2; ++axis) + { + ExpectParticleDomainProjection<2>(axis, isUpper); + } + for (size_t axis = 0; axis < 3; ++axis) + { + ExpectParticleDomainProjection<3>(axis, isUpper); + } + } +} + +TEST(MPMFluidSolver, ParticleColliderProjection) +{ + ExpectParticleColliderProjection<2>(); + ExpectParticleColliderProjection<3>(); +} + +TEST(MPMFluidSolver, HydrostaticColumnRemainsBounded) +{ + ExpectHydrostaticColumnRemainsBounded<2>(); + ExpectHydrostaticColumnRemainsBounded<3>(); +} + +TEST(MPMFluidSolver, DamBreakConservesMassAndStaysInDomain) +{ + ExpectDamBreakConservesMassAndStaysInDomain<2>(); + ExpectDamBreakConservesMassAndStaysInDomain<3>(); +}