From dfd5da287bf261874f08e81d1d7753d950a1e11e Mon Sep 17 00:00:00 2001 From: zhouhang95 <765229842@qq.com> Date: Wed, 24 Nov 2021 23:04:49 +0800 Subject: [PATCH 1/2] axis-bugfix --- zenvis/main.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/zenvis/main.cpp b/zenvis/main.cpp index e060dc3bb0..89e59ca580 100644 --- a/zenvis/main.cpp +++ b/zenvis/main.cpp @@ -68,6 +68,7 @@ void look_perspective( }; grid_blend = ratio_clamp(level - std::floor(level), 0.8, 1.0); center = glm::vec3(0, 0, 0); + radius = 5.0; gizmo_view = glm::lookAt(center - back, center, up); gizmo_proj = glm::ortho(-radius * nx / ny, radius * nx / ny, -radius, radius, -100.0, 100.0); From 600029743aed3256bd40ffa9545dda7ae8a499a6 Mon Sep 17 00:00:00 2001 From: littlemine Date: Fri, 26 Nov 2021 01:21:10 +0800 Subject: [PATCH 2/2] WIP: rebuild mpm pipeline --- projects/gmpm/mpm/Generation.cu | 286 +++++++++++++++++++++++ projects/gmpm/mpm/MPMPipeline.cu | 382 +++++++++++++++++++++++++++++++ projects/gmpm/mpm/MPMStepping.cu | 331 ++++++++++++++++++++++++++ projects/gmpm/mpm/Structures.hpp | 113 +++++++++ 4 files changed, 1112 insertions(+) create mode 100644 projects/gmpm/mpm/Generation.cu create mode 100644 projects/gmpm/mpm/MPMPipeline.cu create mode 100644 projects/gmpm/mpm/MPMStepping.cu create mode 100644 projects/gmpm/mpm/Structures.hpp diff --git a/projects/gmpm/mpm/Generation.cu b/projects/gmpm/mpm/Generation.cu new file mode 100644 index 0000000000..a47c50c89b --- /dev/null +++ b/projects/gmpm/mpm/Generation.cu @@ -0,0 +1,286 @@ +#include "Structures.hpp" + +#include "zensim/cuda/execution/ExecutionPolicy.cuh" +#include "zensim/geometry/VdbSampler.h" +#include "zensim/io/ParticleIO.hpp" +#include "zensim/omp/execution/ExecutionPolicy.hpp" +#include "zensim/tpls/fmt/color.h" +#include "zensim/tpls/fmt/format.h" +#include +#include + +namespace zeno { + +struct ConfigConstitutiveModel : INode { + void apply() override { + auto out = std::make_shared(); + + float dx = get_param("dx"); + if (has_input("dx")) + dx = get_input("dx")->get(); + + float ppc = get_param("ppc"); + if (has_input("ppc")) + ppc = get_input("ppc")->get(); + // volume + out->volume = dx * dx * dx / ppc; + + float density = get_param("density"); + if (has_input("density")) + density = get_input("density")->get(); + // density + out->density = density; + + float E = get_param("E"); + if (has_input("E")) + E = get_input("E")->get(); + + float nu = get_param("nu"); + if (has_input("nu")) + nu = get_input("nu")->get(); + + auto typeStr = get_param("type"); + // elastic model + auto &model = out->getElasticModel(); + if (typeStr == "fcr") + model = zs::FixedCorotated{E, nu}; + else if (typeStr == "nhk") + model = zs::NeoHookean{E, nu}; + else if (typeStr == "stvk") + model = zs::StvkWithHencky{E, nu}; + + // plastic model + out->getPlasticModel() = std::monostate{}; + + set_output("ZSModel", out); + } +}; + +ZENDEFNODE(ConfigConstitutiveModel, { + {"dx", "ppc", "density", "E", "nu"}, + {"ZSModel"}, + {{"float", "dx", "0.1"}, + {"float", "ppc", "8"}, + {"float", "density", "1000"}, + {"string", "type", "fcr"}, + {"float", "E", "10000"}, + {"float", "nu", "0.4"}}, + {"MPM"}, + }); + +struct ToZSParticles : INode { + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing ToZensimParticles\n"); + auto model = get_input("ZSModel"); + + auto inParticles = get_input("prim"); + // const auto size = inParticles->size(); + +#if 0 + auto &obj = inParticles->attr("pos"); +#else + auto obj = zs::sample_from_vdb_file("/home/mine/Codes/Mn/pig.vdb", 0.1, 8); +#endif + + const auto size = obj.size(); + + for (int i = 0; i < obj.size(); ++i) { + auto x = obj[i]; + fmt::print("{}: ({}, {}, {})\n", i, x[0], x[1], x[2]); + if (i % 1000 == 0) + getchar(); + } + + auto outParticles = IObject::make(); + // + outParticles->getModel() = *model; + + // + auto &pars = outParticles->getParticles(); // tilevector + + std::vector tags{ + {"mass", 1}, {"pos", 3}, {"vel", 3}, {"C", 9}, {"vms", 1}}; + + const bool hasPlasticity = model->hasPlasticity(); + const bool hasF = model->hasF(); + + if (hasF) + tags.emplace_back(zs::PropertyTag{"F", 9}); + else + tags.emplace_back(zs::PropertyTag{"J", 1}); + + if (hasPlasticity) + tags.emplace_back(zs::PropertyTag{"logJp", 1}); + + fmt::print("{} particles of these tags\n", size); + for (auto tag : tags) + fmt::print("tag: [{}, {}]\n", tag.name, tag.numChannels); + + { + using namespace zs; + pars = typename ZenoParticles::particles_t{tags, size, memsrc_e::host}; + + auto ompExec = zs::omp_exec(); + ompExec(zs::range(size), [pars = proxy({}, pars), + hasPlasticity, hasF, &inParticles, &model, + &obj](size_t pi) mutable { + using vec3 = zs::vec; + using mat3 = zs::vec; + pars("mass", pi) = model->volume * model->density; + pars.tuple<3>("pos", pi) = obj[pi]; + pars.tuple<3>("vel", pi) = + vec3{0, -1, 0}; // inParticles->attr("vel")[pi]; + pars.tuple<9>("C", pi) = mat3::zeros(); + if (hasF) + pars.tuple<9>("F", pi) = mat3::identity(); + else + pars("J", pi) = 1.; + + if (hasPlasticity) + pars("logJp", pi) = 0; + pars("vms", pi) = 0; // vms + }); + + pars = pars.clone({memsrc_e::um, 0}); + } + + fmt::print(fg(fmt::color::cyan), "done executing ToZensimParticles\n"); + set_output("ZSParticles", outParticles); + } +}; + +ZENDEFNODE(ToZSParticles, { + {"ZSModel", "prim"}, + {"ZSParticles"}, + {}, + {"MPM"}, + }); + +struct MakeZSPartition : INode { + void apply() override { + auto partition = IObject::make(); + partition->get() = + typename ZenoPartition::table_t{(std::size_t)1, zs::memsrc_e::um, 0}; + set_output("ZSPartition", partition); + } +}; +ZENDEFNODE(MakeZSPartition, { + {}, + {"ZSPartition"}, + {}, + {"MPM"}, + }); + +struct MakeZSGrid : INode { + void apply() override { + auto dx = get_param("dx"); + if (has_input("dx")) + dx = get_input("dx")->get(); + + auto grid = IObject::make(); + grid->get() = typename ZenoGrid::grid_t{ + {{"m", 1}, {"v", 3}}, dx, 1, zs::memsrc_e::um, 0}; + + using traits = zs::grid_traits; + fmt::print("grid of dx [{}], side_length [{}], block_size [{}]\n", + grid->get().dx, traits::side_length, traits::block_size); + set_output("ZSGrid", grid); + } +}; +ZENDEFNODE(MakeZSGrid, { + {"dx"}, + {"ZSGrid"}, + {{"float", "dx", "0.1"}}, + {"MPM"}, + }); + +/// conversion + +struct ZSParticlesToPrimitiveObject : zeno::INode { + void apply() override { + fmt::print(fg(fmt::color::green), + "begin executing ZSParticlesToPrimitiveObject\n"); + auto &zspars = get_input("ZSParticles")->getParticles(); + const auto size = zspars.size(); + + auto prim = zeno::IObject::make(); + prim->resize(size); + + using namespace zs; + auto cudaExec = cuda_exec().device(0); + + static_assert(sizeof(zs::vec) == sizeof(zeno::vec3f), + "zeno::vec3f != zs::vec"); + for (auto &&prop : zspars.getPropertyTags()) { + if (prop.numChannels == 3) { + zs::Vector> dst{size, memsrc_e::device, 0}; + cudaExec(zs::range(size), + [zspars = zs::proxy({}, zspars), + dst = zs::proxy(dst), + name = prop.name] __device__(size_t pi) mutable { + dst[pi] = zspars.pack<3>(name, pi); + }); + copy(zs::mem_device, + prim->add_attr(prop.name.asString()).data(), + dst.data(), sizeof(zeno::vec3f) * size); + } else if (prop.numChannels == 1) { + zs::Vector dst{size, memsrc_e::device, 0}; + cudaExec(zs::range(size), + [zspars = zs::proxy({}, zspars), + dst = zs::proxy(dst), + name = prop.name] __device__(size_t pi) mutable { + dst[pi] = zspars(name, pi); + }); + copy(zs::mem_device, prim->add_attr(prop.name.asString()).data(), + dst.data(), sizeof(float) * size); + } + } + fmt::print(fg(fmt::color::cyan), + "done executing ZSParticlesToPrimitiveObject\n"); + set_output("prim", prim); + } +}; + +ZENDEFNODE(ZSParticlesToPrimitiveObject, { + {"ZSParticles"}, + {"prim"}, + {}, + {"MPM"}, + }); + +struct WriteZSParticles : zeno::INode { + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing WriteZSParticles\n"); + auto &pars = get_input("ZSParticles")->getParticles(); + auto path = get_param("path"); + auto cudaExec = zs::cuda_exec().device(0); + zs::Vector> pos{pars.size(), zs::memsrc_e::um, 0}; + zs::Vector vms{pars.size(), zs::memsrc_e::um, 0}; + cudaExec(zs::range(pars.size()), + [pos = zs::proxy(pos), + vms = zs::proxy(vms), + pars = zs::proxy( + {}, pars)] __device__(size_t pi) mutable { + pos[pi] = pars.pack<3>("pos", pi); + vms[pi] = pars("vmstress", pi); + }); + std::vector> posOut(pars.size()); + std::vector vmsOut(pars.size()); + copy(zs::mem_device, posOut.data(), pos.data(), + sizeof(zeno::vec3f) * pars.size()); + copy(zs::mem_device, vmsOut.data(), vms.data(), + sizeof(float) * pars.size()); + + zs::write_partio_with_stress(path, posOut, vmsOut); + fmt::print(fg(fmt::color::cyan), "done executing WriteZSParticles\n"); + } +}; + +ZENDEFNODE(WriteZSParticles, { + {"ZSParticles"}, + {}, + {{"string", "path", ""}}, + {"MPM"}, + }); + +} // namespace zeno \ No newline at end of file diff --git a/projects/gmpm/mpm/MPMPipeline.cu b/projects/gmpm/mpm/MPMPipeline.cu new file mode 100644 index 0000000000..f4add97cee --- /dev/null +++ b/projects/gmpm/mpm/MPMPipeline.cu @@ -0,0 +1,382 @@ +#include "Structures.hpp" + +#include "zensim/cuda/execution/ExecutionPolicy.cuh" +#include "zensim/io/ParticleIO.hpp" +#include "zensim/omp/execution/ExecutionPolicy.hpp" +#include "zensim/simulation/Utils.hpp" +#include "zensim/tpls/fmt/color.h" +#include "zensim/tpls/fmt/format.h" +#include +#include +#include + +namespace zeno { + +/// sparsity +struct ZSPartitionForZSParticles : INode { + void apply() override { + fmt::print(fg(fmt::color::green), + "begin executing ZSPartitionForZSParticles\n"); + auto table = get_input("ZSPartition"); + auto &partition = table->get(); + auto zsgrid = get_input("ZSGrid"); + auto &grid = zsgrid->get(); + + std::vector parObjPtrs{}; + if (has_input("ZSParticles")) + parObjPtrs.push_back(get_input("ZSParticles").get()); + else if (has_input("ZSParticles")) { + auto &objSharedPtrLists = *get_input("ZSParticles"); + for (auto &&objSharedPtr : objSharedPtrLists.get()) + if (auto ptr = dynamic_cast(objSharedPtr.get()); + ptr != nullptr) + parObjPtrs.push_back(ptr); + } + + using namespace zs; + std::size_t cnt = 0; + for (auto &&parObjPtr : parObjPtrs) + cnt += parObjPtr->getParticles().size(); + if (partition._tableSize * 3 / 2 < partition.evaluateTableSize(cnt) || + partition._tableSize / 2 < cnt) + partition.resize(cuda_exec(), cnt); + + using Partition = typename ZenoPartition::table_t; + auto cudaPol = cuda_exec().device(0); + cudaPol(range(partition._tableSize), + [table = proxy(partition)] __device__( + size_t i) mutable { + table._table.keys[i] = + Partition::key_t::uniform(Partition::key_scalar_sentinel_v); + table._table.indices[i] = Partition::sentinel_v; + table._table.status[i] = -1; + if (i == 0) + *table._cnt = 0; + }); + using grid_t = typename ZenoGrid::grid_t; + static_assert(grid_traits::is_power_of_two, + "grid side_length should be power of two"); + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + cudaPol(range(pars.size()), + [pars = proxy({}, pars), + table = proxy(partition), + dxinv = 1.f / grid.dx] __device__(size_t pi) mutable { + auto x = pars.template pack<3>("pos", pi); + auto c = (x * dxinv - 0.5); + typename Partition::key_t coord{}; + for (int d = 0; d != 3; ++d) + coord[d] = lower_trunc(c[d]); + table.insert(coord - (coord & (grid_t::side_length - 1))); + }); + } + fmt::print("partition of [{}] blocks for {} particles\n", partition.size(), + cnt); + + fmt::print(fg(fmt::color::cyan), + "done executing ZSPartitionForZSParticles\n"); + set_output("ZSPartition", table); + } +}; + +ZENDEFNODE(ZSPartitionForZSParticles, + { + {"ZSPartition", "ZSGrid", "ZSParticles"}, + {"ZSPartition"}, + {}, + {"MPM"}, + }); + +struct ExpandZSPartition : INode { + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing ExpandZSPartition\n"); + auto table = get_input("ZSPartition"); + auto &partition = table->get(); + auto offset = get_param("offset"); + auto extent = get_param("extent"); + + using namespace zs; + auto cudaPol = cuda_exec().device(0); + using grid_t = typename ZenoGrid::grid_t; + static_assert(grid_traits::is_power_of_two, + "grid side_length should be power of two"); + + auto prevCnt = partition.size(); + cudaPol(range(prevCnt), [table = proxy(partition), + offset, extent] __device__(size_t bi) mutable { + auto blockid = table._activeKeys[bi]; + for (auto ijk : ndrange<3>(extent)) + table.insert(blockid + (make_vec(ijk) + offset) * + (int)grid_traits::side_length); + }); + fmt::print("partition insertion [{}] blocks -> [{}] blocks\n", prevCnt, + partition.size()); + fmt::print(fg(fmt::color::cyan), "done executing ExpandZSPartition\n"); + + set_output("ZSPartition", std::move(table)); + } +}; + +ZENDEFNODE(ExpandZSPartition, + { + {"ZSPartition"}, + {"ZSPartition"}, + {{"int", "offset", "0"}, {"int", "extent", "2"}}, + {"MPM"}, + }); + +/// grid +struct ZSGridFromZSPartition : INode { + void apply() override { + fmt::print(fg(fmt::color::green), + "begin executing ZSGridFromZSPartition\n"); + auto &partition = get_input("ZSPartition")->get(); + auto cnt = partition.size(); + + auto zsgrid = get_input("ZSGrid"); + auto &grid = zsgrid->get(); + grid.resize(cnt); + + using namespace zs; + auto cudaPol = cuda_exec().device(0); + // clear grid + cudaPol({(int)cnt, (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid)] __device__( + int bi, int ci) mutable { + auto block = grid.block(bi); + const auto nchns = grid.numChannels(); + for (int i = 0; i != nchns; ++i) + block(i, ci) = 0; + }); + + fmt::print(fg(fmt::color::cyan), "done executing ZSGridFromZSPartition\n"); + set_output("ZSGrid", zsgrid); + } +}; + +ZENDEFNODE(ZSGridFromZSPartition, { + {"ZSPartition", "ZSGrid"}, + {"ZSGrid"}, + {}, + {"MPM"}, + }); + +struct UpdateZSGrid : INode { + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing UpdateZSGrid\n"); + // auto dt = get_input("dt")->as()->get(); + auto maxVelSqr = IObject::make(); + + auto &partition = get_input("ZSPartition")->get(); + auto zsgrid = get_input("ZSGrid"); + auto &grid = zsgrid->get(); + auto stepDt = get_input("dt")->get(); + + using namespace zs; + auto gravity = get_param("gravity"); + auto accel = zs::vec::zeros(); + if (has_input("ExtForce")) { + auto tmp = get_input("Accel")->get(); + accel = zs::vec{tmp[0], tmp[1], tmp[2]}; + } else + accel[1] = gravity; + + Vector velSqr{1, zs::memsrc_e::um, 0}; + velSqr[0] = 0; + auto cudaPol = cuda_exec().device(0); + + cudaPol({(int)partition.size(), (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid), stepDt, accel, + ptr = velSqr.data()] __device__(int bi, int ci) mutable { + auto block = grid.block(bi); + auto mass = block("m", ci); + if (mass != 0.f) { + mass = 1.f / mass; + auto vel = block.pack<3>("v", ci) * mass; + vel += accel * stepDt; + block.set("v", ci, vel); + /// cfl dt + auto velSqr = vel.l2NormSqr(); + atomic_max(exec_cuda, ptr, velSqr); + } + }); + + maxVelSqr->set(velSqr[0]); + fmt::print(fg(fmt::color::cyan), "done executing GridUpdate\n"); + set_output("ZSGrid", zsgrid); + set_output("MaxVelSqr", maxVelSqr); + } +}; + +ZENDEFNODE(UpdateZSGrid, { + {"ZSPartition", "ZSGrid", "dt", "Accel"}, + {"ZSGrid", "MaxVelSqr"}, + {{"float", "gravity", "-9.8"}}, + {"MPM"}, + }); + +struct ZSParticleToZSGrid : INode { + template + void p2g(zs::CudaExecutionPolicy &cudaPol, const Model &model, + const float volume, const typename ZenoParticles::particles_t &pars, + const typename ZenoPartition::table_t &partition, const float dt, + typename ZenoGrid::grid_t &grid) { + using namespace zs; + cudaPol(range(pars.size()), [pars = proxy({}, pars), + table = proxy(partition), + grid = proxy({}, grid), dt, + dxinv = 1.f / grid.dx, vol = volume, + model] __device__(size_t pi) mutable { + using grid_t = RM_CVREF_T(grid); + const auto Dinv = 4.f * dxinv * dxinv; + auto localPos = pars.pack<3>("pos", pi); + auto vel = pars.pack<3>("vel", pi); + auto mass = pars("mass", pi); + auto C = pars.pack<3, 3>("C", pi); + auto F = pars.pack<3, 3>("F", pi); + auto P = model.first_piola(F); + + auto contrib = -dt * Dinv * vol * P * F.transpose(); + auto arena = make_local_arena(grid.dx, localPos); + + for (auto loc : arena.range()) { + auto coord = arena.coord(loc); + auto localIndex = coord & (grid_t::side_length - 1); + auto blockno = table.query(coord - localIndex); + if (blockno < 0) + printf("THE HELL!"); + auto block = grid.block(blockno); + + auto xixp = arena.diff(loc); + auto W = arena.weight(loc); + const auto cellid = grid_t::coord_to_cellid(localIndex); + atomic_add(exec_cuda, &block("m", cellid), mass * W); + auto Cxixp = C * xixp; + auto fdt = contrib * xixp; + for (int d = 0; d != 3; ++d) + atomic_add(exec_cuda, &block("v", d, cellid), + W * (mass * (vel[d] + Cxixp[d]) + fdt[d])); + } + }); + } + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing ZSParticleToZSGrid\n"); + + std::vector parObjPtrs{}; + if (has_input("ZSParticles")) + parObjPtrs.push_back(get_input("ZSParticles").get()); + else if (has_input("ZSParticles")) { + auto &objSharedPtrLists = *get_input("ZSParticles"); + for (auto &&objSharedPtr : objSharedPtrLists.get()) + if (auto ptr = dynamic_cast(objSharedPtr.get()); + ptr != nullptr) + parObjPtrs.push_back(ptr); + } + auto &partition = get_input("ZSPartition")->get(); + auto zsgrid = get_input("ZSGrid"); + auto &grid = zsgrid->get(); + auto stepDt = get_input("dt")->get(); + + using namespace zs; + auto cudaPol = cuda_exec().device(0); + + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + auto &model = parObjPtr->getModel(); + + fmt::print("[p2g] dx: {}, dt: {}, npars: {}\n", grid.dx, stepDt, + pars.size()); + + match([&](auto &elasticModel) { + p2g(cudaPol, elasticModel, model.volume, pars, partition, stepDt, grid); + })(model.getElasticModel()); + } + + fmt::print(fg(fmt::color::cyan), "done executing ZSParticleToZSGrid\n"); + set_output("ZSGrid", zsgrid); + } +}; + +ZENDEFNODE(ZSParticleToZSGrid, + { + {"ZSParticles", "ZSPartition", "ZSGrid", "dt"}, + {"ZSGrid"}, + {}, + {"MPM"}, + }); + +struct ZSGridToZSParticle : INode { + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing ZSGridToZSParticle\n"); + auto &grid = get_input("ZSGrid")->get(); + auto &partition = get_input("ZSPartition")->get(); + + std::vector parObjPtrs{}; + if (has_input("ZSParticles")) + parObjPtrs.push_back(get_input("ZSParticles").get()); + else if (has_input("ZSParticles")) { + auto &objSharedPtrLists = *get_input("ZSParticles"); + for (auto &&objSharedPtr : objSharedPtrLists.get()) + if (auto ptr = dynamic_cast(objSharedPtr.get()); + ptr != nullptr) + parObjPtrs.push_back(ptr); + } + + auto stepDt = get_input("dt")->get(); + + using namespace zs; + auto cudaPol = cuda_exec().device(0); + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + + cudaPol(range(pars.size()), + [pars = proxy({}, pars), + table = proxy(partition), + grid = proxy({}, grid), dt = stepDt, + dxinv = 1.f / grid.dx] __device__(size_t pi) mutable { + using grid_t = RM_CVREF_T(grid); + const auto Dinv = 4.f * dxinv * dxinv; + auto pos = pars.pack<3>("pos", pi); + auto vel = zs::vec::zeros(); + auto C = zs::vec::zeros(); + + auto arena = make_local_arena(grid.dx, pos); + for (auto loc : arena.range()) { + auto coord = arena.coord(loc); + auto localIndex = coord & (grid_t::side_length - 1); + auto blockno = table.query(coord - localIndex); + if (blockno < 0) + printf("THE HELL!"); + auto block = grid.block(blockno); + auto xixp = arena.diff(loc); + auto W = arena.weight(loc); + auto vi = + block.pack<3>("v", grid_t::coord_to_cellid(localIndex)); + + vel += vi * W; + C += W * Dinv * dyadic_prod(vi, xixp); + } + pars.tuple<3>("vel", pi) = vel; + pars.tuple<3 * 3>("C", pi) = C; + pos += vel * dt; + pars.tuple<3>("pos", pi) = pos; + + auto F = pars.pack<3, 3>("F", pi); + auto tmp = zs::vec::identity() + C * dt; + F = tmp * F; + pars.tuple<3 * 3>("F", pi) = F; + }); + } + fmt::print(fg(fmt::color::cyan), "done executing ZSGridToZSParticle\n"); + } +}; + +ZENDEFNODE(ZSGridToZSParticle, + { + {"ZSGrid", "ZSPartition", "ZSParticles", "dt"}, + {}, + {}, + {"MPM"}, + }); + +} // namespace zeno \ No newline at end of file diff --git a/projects/gmpm/mpm/MPMStepping.cu b/projects/gmpm/mpm/MPMStepping.cu new file mode 100644 index 0000000000..d063618837 --- /dev/null +++ b/projects/gmpm/mpm/MPMStepping.cu @@ -0,0 +1,331 @@ +#include "Structures.hpp" + +#include "zensim/cuda/execution/ExecutionPolicy.cuh" +#include "zensim/io/ParticleIO.hpp" +#include "zensim/omp/execution/ExecutionPolicy.hpp" +#include "zensim/simulation/Utils.hpp" +#include "zensim/tpls/fmt/color.h" +#include "zensim/tpls/fmt/format.h" +#include +#include +#include + +namespace zeno { + +struct MPMStepping : INode { + template + void p2g(zs::CudaExecutionPolicy &cudaPol, const Model &model, + const float volume, const typename ZenoParticles::particles_t &pars, + const typename ZenoPartition::table_t &partition, const float dt, + typename ZenoGrid::grid_t &grid) { + using namespace zs; + cudaPol(range(pars.size()), [pars = proxy({}, pars), + table = proxy(partition), + grid = proxy({}, grid), dt, + dxinv = 1.f / grid.dx, vol = volume, + model] __device__(size_t pi) mutable { + using grid_t = RM_CVREF_T(grid); + const auto Dinv = 4.f * dxinv * dxinv; + auto localPos = pars.pack<3>("pos", pi); + auto vel = pars.pack<3>("vel", pi); + auto mass = pars("mass", pi); + auto C = pars.pack<3, 3>("C", pi); + auto F = pars.pack<3, 3>("F", pi); + auto P = model.first_piola(F); + + auto contrib = -dt * Dinv * vol * P * F.transpose(); + auto arena = make_local_arena(grid.dx, localPos); + + if (pars.size() < 3000 && pi < 3) { + auto p = localPos; + auto v = vel; + printf("[p2g] pos[%d] (%f, %f, %f) vel (%f, %f, %f) mass (%f)\n", + (int)pi, p[0], p[1], p[2], v[0], v[1], v[2], mass); + } + + for (auto loc : arena.range()) { + auto coord = arena.coord(loc); + auto localIndex = coord & (grid_t::side_length - 1); + auto blockno = table.query(coord - localIndex); + if (blockno < 0) + printf("THE HELL!"); + auto block = grid.block(blockno); + + auto xixp = arena.diff(loc); + auto W = arena.weight(loc); + const auto cellid = grid_t::coord_to_cellid(localIndex); + atomic_add(exec_cuda, &block("m", cellid), mass * W); + auto Cxixp = C * xixp; + auto fdt = contrib * xixp; + for (int d = 0; d != 3; ++d) + atomic_add(exec_cuda, &block("v", d, cellid), + W * (mass * (vel[d] + Cxixp[d]) + fdt[d])); + } + }); + } + void apply() override { + fmt::print(fg(fmt::color::green), "begin executing MPMStepping\n"); + + size_t numPars = 0; + std::vector parObjPtrs{}; + if (has_input("ZSParticles")) + parObjPtrs.push_back(get_input("ZSParticles").get()); + else if (has_input("ZSParticles")) { + auto &objSharedPtrLists = *get_input("ZSParticles"); + for (auto &&objSharedPtr : objSharedPtrLists.get()) + if (auto ptr = dynamic_cast(objSharedPtr.get()); + ptr != nullptr) + parObjPtrs.push_back(ptr); + } + + for (auto parObjPtr : parObjPtrs) + numPars += parObjPtr->getParticles().size(); + + auto &partition = get_input("ZSPartition")->get(); + auto zsgrid = get_input("ZSGrid"); + auto &grid = zsgrid->get(); + auto stepDt = get_input("dt")->get(); + + using namespace zs; + auto cudaPol = cuda_exec().device(0); + +#if 0 + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + cudaPol(range(pars.size()), + [pars = proxy({}, pars), + table = proxy(partition), + grid = proxy({}, grid), dt = stepDt, + dxinv = 1.f / grid.dx] __device__(size_t pi) mutable { + using grid_t = RM_CVREF_T(grid); + auto pos = pars.pack<3>("pos", pi); + auto vel = pars.pack<3>("vel", pi); + pos += vel * dt; + pars.tuple<3>("pos", pi) = pos; + }); + } +#endif + + /// partition + using Partition = typename ZenoPartition::table_t; + using grid_t = typename ZenoGrid::grid_t; + partition.resize(cudaPol, numPars); + cudaPol(range(partition._tableSize), + [table = proxy(partition)] __device__( + size_t i) mutable { + table._table.keys[i] = + Partition::key_t::uniform(Partition::key_scalar_sentinel_v); + table._table.indices[i] = Partition::sentinel_v; + table._table.status[i] = -1; + if (i == 0) + *table._cnt = 0; + }); + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + cudaPol(range(pars.size()), + [pars = proxy({}, pars), + table = proxy(partition), + dxinv = 1.f / grid.dx] __device__(size_t pi) mutable { + auto x = pars.pack<3>(1, pi); + auto c = (x * dxinv - 0.5); + typename Partition::key_t coord{}; + for (int d = 0; d != 3; ++d) + coord[d] = lower_trunc(c[d]); + table.insert((coord - (coord & (grid_t::side_length - 1)))); + }); + } + cudaPol(range(partition.size()), + [table = proxy(partition), offset = 0, + extent = 2] __device__(size_t bi) mutable { + auto blockid = table._activeKeys[bi]; + for (auto ijk : ndrange<3>(extent)) + table.insert(blockid + (make_vec(ijk) + offset) * + grid_t::side_length); + }); + /// grid + grid.resize(1000000); + // clear grid + fmt::print("blocksize: {}\n", ZenoGrid::grid_t::block_space()); + cudaPol({(int)partition.size(), (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid)] __device__( + int bi, int ci) mutable { + auto block = grid.block(bi); + const auto nchns = grid.numChannels(); + for (int i = 0; i != nchns; ++i) + block(i, ci) = 0; + }); + + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + auto &model = parObjPtr->getModel(); + + fmt::print("[p2g] dx: {}, dt: {}, npars: {}\n", grid.dx, stepDt, + pars.size()); + + match([&](auto &elasticModel) { + p2g(cudaPol, elasticModel, model.volume, pars, partition, stepDt, grid); + })(model.getElasticModel()); + } + + { + /// particle momentum + Vector lm{3, memsrc_e::um, 0}, am{3, memsrc_e::um, 0}; + for (int d = 0; d != 3; ++d) + lm[d] = am[d] = 0; + auto &pars = parObjPtrs[0]->getParticles(); + cudaPol(range(pars.size()), [pars = proxy({}, pars), + lm = proxy(lm), + am = proxy(am), + dx = grid.dx] __device__(size_t pi) { + auto mass = pars("mass", pi); + auto pos = pars.pack<3>("pos", pi); + auto vel = pars.pack<3>("vel", pi); + const auto Dinv = 4.f / dx / dx; + auto B = pars.pack<3, 3>("C", pi) / Dinv; + for (int d = 0; d != 3; ++d) + atomic_add(exec_cuda, &lm[d], mass * vel[d]); + + auto res = pos.cross(vel) * mass; + for (int a = 0; a != 3; ++a) + for (int b = 0; b != 3; ++b) + for (int g = 0; g != 3; ++g) + if ((a == 0 && b == 1 && g == 2) || + (a == 1 && b == 2 && g == 0) || (a == 2 && b == 0 && g == 1)) + res(g) += mass * B(a, b); + else if ((a == 0 && b == 2 && g == 1) || + (a == 1 && b == 0 && g == 2) || + (a == 2 && b == 1 && g == 0)) + res(g) -= mass * B(a, b); + + for (int d = 0; d != 3; ++d) + atomic_add(exec_cuda, &am[d], res[d]); + }); + fmt::print("[[after p2g]] pars({} total) linear momentum: ({}, {}, {}), " + "angular momentum: ({}, {}, {})\n", + pars.size(), lm[0], lm[1], lm[2], am[0], am[1], am[2]); + } + + auto cnt = partition.size(); + cudaPol({(int)partition.size(), (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid), + stepDt] __device__(int bi, int ci) mutable { + auto block = grid.block(bi); + auto mass = block("m", ci); + if (mass != 0.f) { + mass = 1.f / mass; + auto vel = block.pack<3>("v", ci) * mass; + // vel += extf * stepDt; + // vel[1] += -9.8 * stepDt; + /// write back + block.set("v", ci, vel); + } + }); + fmt::print("grid update numblocks {}\n", cnt); + + { + Vector lm{3, memsrc_e::um, 0}, am{3, memsrc_e::um, 0}; + /// particle momentum + for (int d = 0; d != 3; ++d) + lm[d] = am[d] = 0; + cudaPol({(int)partition.size(), (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid), + table = proxy(partition), + lm = proxy(lm), + am = proxy(am)] __device__(int bi, int ci) { + using grid_t = RM_CVREF_T(grid); + auto block = grid.block(bi); + auto mass = block("m", ci); + if (mass != 0) { + auto blockkey = table._activeKeys[bi]; + auto x = (blockkey + grid_t::cellid_to_coord(ci)) * grid.dx; + auto mv = block.pack<3>("v", ci) * mass; + /// x cross mv; + auto res = x.cross(mv); + + for (int i = 0; i != 3; ++i) { + atomic_add(exec_cuda, &am[i], res[i]); + atomic_add(exec_cuda, &lm[i], mv[i]); + } + } + }); + fmt::print("[[after grid update]] grid({} blocks total) linear momentum: " + "({}, {}, {}), angular " + "momentum: ({}, {}, {})\n", + partition.size(), lm[0], lm[1], lm[2], am[0], am[1], am[2]); + } + + Vector vols{2, memsrc_e::um, 0}; + vols[0] = vols[1] = 0.; + cudaPol({(int)partition.size(), (int)ZenoGrid::grid_t::block_space()}, + [grid = proxy({}, grid), + vols = proxy(vols), + stepDt] __device__(int bi, int ci) mutable { + auto block = grid.block(bi); + auto mass = block("m", ci); + if (mass != 0.f) { + atomic_add(exec_cuda, &vols[0], + (double)(grid.dx * grid.dx * grid.dx)); + atomic_add(exec_cuda, &vols[1], 1.); + } + }); + + fmt::print(fg(fmt::color::red), + "sum total particle volume: {} x {} => {}\n", numPars, + parObjPtrs[0]->getModel().volume, + parObjPtrs[0]->getModel().volume * numPars); + fmt::print(fg(fmt::color::brown), "sum total grid volume: {} x {} => {}\n", + vols[1], grid.dx * grid.dx * grid.dx, vols[0]); + + for (auto &&parObjPtr : parObjPtrs) { + auto &pars = parObjPtr->getParticles(); + cudaPol(range(pars.size()), + [pars = proxy({}, pars), + table = proxy(partition), + grid = proxy({}, grid), dt = stepDt, + dxinv = 1.f / grid.dx] __device__(size_t pi) mutable { + using grid_t = RM_CVREF_T(grid); + const auto Dinv = 4.f * dxinv * dxinv; + auto pos = pars.pack<3>("pos", pi); + auto vel = zs::vec::zeros(); + auto C = zs::vec::zeros(); + + auto arena = make_local_arena(grid.dx, pos); + for (auto loc : arena.range()) { + auto coord = arena.coord(loc); + auto localIndex = coord & (grid_t::side_length - 1); + auto blockno = table.query(coord - localIndex); + if (blockno < 0) + printf("THE HELL!"); + auto block = grid.block(blockno); + auto xixp = arena.diff(loc); + auto W = arena.weight(loc); + auto vi = + block.pack<3>("v", grid_t::coord_to_cellid(localIndex)); + + vel += vi * W; + C += W * Dinv * dyadic_prod(vi, xixp); + } + pars.tuple<3>("vel", pi) = vel; + pars.tuple<3 * 3>("C", pi) = C; + pos += vel * dt; + pars.tuple<3>("pos", pi) = pos; + + auto F = pars.pack<3, 3>("F", pi); + auto tmp = zs::vec::identity() + C * dt; + F = tmp * F; + pars.tuple<3 * 3>("F", pi) = F; + }); + } + + fmt::print(fg(fmt::color::cyan), "done executing MPMStepping\n"); + } +}; + +ZENDEFNODE(MPMStepping, { + {"ZSParticles", "ZSPartition", "ZSGrid", "dt"}, + {}, + {}, + {"MPM"}, + }); + +} // namespace zeno \ No newline at end of file diff --git a/projects/gmpm/mpm/Structures.hpp b/projects/gmpm/mpm/Structures.hpp new file mode 100644 index 0000000000..9413527d68 --- /dev/null +++ b/projects/gmpm/mpm/Structures.hpp @@ -0,0 +1,113 @@ +#pragma once +#include "zensim/container/HashTable.hpp" +#include "zensim/container/IndexBuckets.hpp" +#include "zensim/geometry/AnalyticLevelSet.h" +#include "zensim/geometry/Collider.h" +#include "zensim/geometry/SparseLevelSet.hpp" +#include "zensim/geometry/Structure.hpp" +#include "zensim/geometry/Structurefree.hpp" +#include "zensim/physics/ConstitutiveModel.hpp" +#include "zensim/physics/constitutive_models/EquationOfState.hpp" +#include "zensim/physics/constitutive_models/FixedCorotated.h" +#include "zensim/physics/constitutive_models/NeoHookean.hpp" +#include "zensim/physics/constitutive_models/StvkWithHencky.hpp" +#include "zensim/physics/plasticity_models/VonMisesCapped.hpp" +#include "zensim/resource/Resource.h" +#include + +namespace zeno { + +using ElasticModel = + zs::variant, zs::NeoHookean, + zs::StvkWithHencky>; +using PlasticModel = zs::variant; + +struct ZenoConstitutiveModel : IObject { + enum elastic_model_e { Fcr, Nhk, Stvk }; + enum plastic_model_e { None, VonMises, DruckerPrager, Camclay }; + + auto &getElasticModel() noexcept { return elasticModel; } + const auto &getElasticModel() const noexcept { return elasticModel; } + auto &getPlasticModel() noexcept { return plasticModel; } + const auto &getPlasticModel() const noexcept { return plasticModel; } + + template auto &getElasticModel() noexcept { + return std::get(elasticModel); + } + template const auto &getElasticModel() const noexcept { + return std::get(elasticModel); + } + template auto &getPlasticModel() noexcept { + return std::get(plasticModel); + } + template const auto &getPlasticModel() const noexcept { + return std::get(plasticModel); + } + + bool hasF() const noexcept { return elasticModel.index() < 3; } + bool hasPlasticity() const noexcept { return plasticModel.index() != 0; } + + float volume, density; + ElasticModel elasticModel; + PlasticModel plasticModel; +}; + +struct ZenoParticles : IObject { + using particles_t = + zs::TileVector>; + auto &getParticles() noexcept { return particles; } + const auto &getParticles() const noexcept { return particles; } + auto &getModel() noexcept { return model; } + const auto &getModel() const noexcept { return model; } + particles_t particles{}; + ZenoConstitutiveModel model{}; +}; + +struct ZenoGrid : IObject { + using grid_t = + zs::Grid>; + auto &get() noexcept { return grid; } + const auto &get() const noexcept { return grid; } + grid_t grid; +}; + +struct ZenoPartition : IObject { + using table_t = zs::HashTable>; + auto &get() noexcept { return table; } + const auto &get() const noexcept { return table; } + table_t table; +}; + +struct ZenoIndexBuckets : IObject { + using buckets_t = zs::IndexBuckets<3, int, int>; + auto &get() noexcept { return ibs; } + const auto &get() const noexcept { return ibs; } + buckets_t ibs; +}; + +struct ZenoSparseLevelSet : IObject { + using spls_t = zs::SparseLevelSet<3>; + auto &get() noexcept { return ls; } + const auto &get() const noexcept { return ls; } + spls_t ls; +}; + +struct ZenoBoundary : IObject { + using boundary_t = + zs::variant>, + zs::Collider>>; + enum category_e { LevelSet, Plane }; + + auto &getBoundary() noexcept { return boundary; } + const auto &getBoundary() const noexcept { return boundary; } + template auto &getBoundary() noexcept { + return std::get(boundary); + } + template const auto &getBoundary() const noexcept { + return std::get(boundary); + } + boundary_t boundary; +}; + +} // namespace zeno \ No newline at end of file