From 7608b7a13ddc16b9b9838e50166dec89dfb3afac Mon Sep 17 00:00:00 2001 From: cflin Date: Tue, 11 Apr 2023 21:37:12 +0800 Subject: [PATCH] add boundary visible && mesh visible options --- CMakeLists.txt | 1 + src/viewer/UIStaticMenu.cpp | 31 +++-- src/viewer/UIStaticSimState.h | 103 ++++++++++++++-- static_sim/StaticSim.cpp | 20 +-- static_sim/StaticSim.h | 224 ++++++++++++++++++++-------------- 5 files changed, 254 insertions(+), 125 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2cedeb5..0f6f09a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -367,3 +367,4 @@ add_subdirectory(sim-test/rigid-test) # IPC link static simulation target_link_libraries(rigid_ipc_sim PUBLIC StaticSim) +add_compile_definitions(DEBUG_SSIM) \ No newline at end of file diff --git a/src/viewer/UIStaticMenu.cpp b/src/viewer/UIStaticMenu.cpp index 2095cac..8008b0e 100644 --- a/src/viewer/UIStaticMenu.cpp +++ b/src/viewer/UIStaticMenu.cpp @@ -83,10 +83,14 @@ namespace ipc::rigid { void UIStaticSimState::draw_io() { if (ImGui::Button("加载文件##IO", ImVec2(-1, 0))) { +#ifndef DEBUG_SSIM std::string fname = igl::file_dialog_open(); if (fname != "") { load(fname); } +#else + load("/home/cflin/Documents/CppField/source/RigidElasticSim/sim-test/rigid-test/cube/config.json"); +#endif } if (m_has_scene) { // if (ImGui::Button("重新加载##IO", ImVec2(-1, 0))) { @@ -149,21 +153,28 @@ namespace ipc::rigid { } // ImGui::RadioButton("开始计算##SimPlayer", &player_state, PlayerState::Playing); // ImGui::RadioButton("暂停##SimPlayer", &player_state, PlayerState::Paused); - + ImGui::Checkbox("边界条件", &gui_ctrl_.is_visible_BC); + ImGui::Checkbox("网格可视化",&gui_ctrl_.mesh_visible); if (ImGui::CollapsingHeader( "边界条件设置##SimPlayer", ImGuiTreeNodeFlags_DefaultOpen)) { - // -------------------------------------------------------------------- ImGui::Checkbox("边界包围盒", &m_bkp_has_intersections); + // -------------------------------------------------------------------- - // ImGui::Checkbox("受力", &m_state.m_solve_collisions); - static float zoom = 50.0f; - static float zoom1 = 50.0f; - static float zoom2 = 50.0f; - ImGui::SliderFloat("x-受力", &zoom, -50, 50, "%.3f"); - ImGui::SliderFloat("y-受力", &zoom1, -50, 50, "%.3f"); - ImGui::SliderFloat("z-受力", &zoom2, -50, 50, "%.3f"); - } + for (int i = 0; i < sp_StaticSim_->DirichletBCs.size(); ++i) { + // each boundary + // ImGui::Checkbox("受力", &m_state.m_solve_collisions); + auto &i_DBC = sp_StaticSim_->DirichletBCs[i]; + std::vector v_min_point = {(float) i_DBC.relMinBBox.x(), (float) i_DBC.relMinBBox.y(), + (float) i_DBC.relMinBBox.z()}; + std::vector v_max_point = {(float) i_DBC.relMaxBBox.x(), (float) i_DBC.relMaxBBox.y(), + (float) i_DBC.relMaxBBox.z()}; + ImGui::SliderFloat3("最小值", v_min_point.data(), -50, 50, "%.3f"); + ImGui::SliderFloat3("最大值", v_max_point.data(), -50, 50, "%.3f"); + i_DBC.relMinBBox={v_min_point[0],v_min_point[1],v_min_point[2]}; + i_DBC.relMaxBBox={v_max_point[0],v_max_point[1],v_max_point[2]}; + } + } // ImGui::Checkbox("受力大小", &m_bkp_optimization_failed); // ImGui::SameLine(); // ImGui::HelpMarker("yes - stop playing if step optimization failed."); diff --git a/src/viewer/UIStaticSimState.h b/src/viewer/UIStaticSimState.h index 4e71c54..812170c 100644 --- a/src/viewer/UIStaticSimState.h +++ b/src/viewer/UIStaticSimState.h @@ -16,9 +16,7 @@ #include #include #include -#include -#include -#include +//#include #include "UISimState.hpp" @@ -33,6 +31,7 @@ namespace ipc { bool is_modified = false; bool is_solved = false; bool is_loaded_json = false; + bool is_visible_BC = false; bool mesh_visible = true; bool single_color = false; ssim::SimTargetOption::Target target_to_evaluate = ssim::SimTargetOption::U_NORM; @@ -63,6 +62,12 @@ namespace ipc { // back_color m_viewer.core().background_color << 0.9f, 0.9f, 0.9f, 0.4f; + for (int i = 0; i < 1 + sp_StaticSim_->DirichletBCs.size() + sp_StaticSim_->NeumannBCs.size(); ++i) { + viewer->append_mesh(); + if (i > 0) { + viewer->data(i).show_lines = false; + } + } } @@ -75,7 +80,8 @@ namespace ipc { const Eigen::MatrixXi &F = model.F; // plot - viewer->data().set_mesh(V, F); + viewer->data(0).set_mesh(V, F); +// viewer->data().add_edges(V.topRows(10),V.topRows(10).array()+10,Eigen::RowVector3d(0,0,1)); // viewer->data().set_colors(Eigen::Vector3d(0.5,0.5,1)); m_has_scene = true; @@ -131,6 +137,7 @@ namespace ipc { // void load_scene(); void redraw_scene() override { + // get V, F ssim::Model model = sp_StaticSim_->get_mesh(); const Eigen::MatrixXd &V = model.V; const Eigen::MatrixXi &F = model.F; @@ -138,32 +145,105 @@ namespace ipc { Eigen::MatrixXd target = sp_StaticSim_->EvaluateTarget(gui_ctrl_.target_to_evaluate); // plot - viewer->data().set_mesh(V, F); + viewer->data(0).set_mesh(V, F); + + + viewer->data(0).show_lines = gui_ctrl_.mesh_visible; - if (!gui_ctrl_.mesh_visible) { - viewer->data().show_overlay = false; - viewer->data().show_lines = false; - } if (!gui_ctrl_.single_color) { // colormap Eigen::MatrixXd C; igl::jet(target, true, C); - viewer->data().set_colors(C); + viewer->data(0).set_colors(C); } else { // single color - viewer->data().set_colors(Eigen::RowVector3d(1, 0, 0)); + viewer->data(0).set_colors(Eigen::RowVector3d(1, 0, 0)); // colorbar_plugin_.draw_colorbar(igl::ColorMapType::COLOR_MAP_TYPE_JET,0,1,/*background_color=*/Eigen::Vector4f(1,1,1,0)); } } + void load_bc() { + auto get_V_F_from_bd_box = [](const Eigen::Vector3d &v_min, + const Eigen::Vector3d &v_max) -> std::tuple { + Eigen::MatrixXd V(8, 3); + V << -1, -1, -1, + 1, -1, -1, + 1, 1, -1, + -1, 1, -1, + -1, -1, 1, + 1, -1, 1, + 1, 1, 1, + -1, 1, 1; + // Define the faces of the cube + static Eigen::MatrixXi F(12, 3); + F << 0, 2, 1, + 0, 3, 2, + 4, 5, 6, + 4, 6, 7, + 3, 6, 2, + 3, 7, 6, + 0, 1, 5, + 0, 5, 4, + 5, 1, 2, + 5, 2, 6, + 4, 7, 3, + 5, 4, 3; + Eigen::Vector3d transform = (v_max + v_min) / 2.0; + Eigen::Vector3d scaled_factor = (v_max - v_min) / 2.0; + V = V.array().rowwise() * scaled_factor.transpose().array(); + V = V.array().rowwise() + transform.transpose().array(); + return {V, F}; + + }; + +// // get model size +// Eigen::Vector3d model_size = sp_StaticSim_->get_mesh().V.colwise().maxCoeff() - +// sp_StaticSim_->get_mesh().V.colwise().minCoeff(); +// double longest_size = model_size.maxCoeff(); +//// spdlog::warn("longest_size: {}", longest_size); +// +// // DBC, NBC +// ssim::Model mesh_in_DBC = sp_StaticSim_->GetMeshInDBC(); +// // Normal +// Eigen::MatrixXd N; +// igl::per_vertex_normals(mesh_in_DBC.V, mesh_in_DBC.F, N); +//// spdlog::warn("VFN: {}, {}, {}", mesh_in_DBC.V.rows(), mesh_in_DBC.F.rows(), N.rows()); +// +// +// // Draw arrows toward the vertex normals +// Eigen::MatrixXd starts = mesh_in_DBC.V; +// Eigen::MatrixXd ends = starts.array() + 0.1 * longest_size * N.array(); +// viewer->data().add_edges(starts, ends, Eigen::RowVector3d(0, 0, 1)); + for (int i = 0; i < sp_StaticSim_->DirichletBCs.size(); ++i) { + auto [V, F] = get_V_F_from_bd_box(sp_StaticSim_->DirichletBCs[i].relMinBBox, + sp_StaticSim_->DirichletBCs[i].relMaxBBox); + viewer->data(i + 1).set_mesh(V, F); + viewer->data(i + 1).set_colors(Eigen::RowVector3d(1, 0, 0)); + viewer->data(i + 1).is_visible = gui_ctrl_.is_visible_BC; + + } + for (int i = 0; i < sp_StaticSim_->NeumannBCs.size(); ++i) { + auto [V, F] = get_V_F_from_bd_box(sp_StaticSim_->NeumannBCs[i].relMinBBox, + sp_StaticSim_->NeumannBCs[i].relMaxBBox); + viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).set_mesh(V, F); + viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).set_colors(Eigen::RowVector3d(0, 0, 1)); + viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).is_visible = gui_ctrl_.is_visible_BC; + + } + + } + bool pre_draw_loop() { if (gui_ctrl_.is_loaded_json && !gui_ctrl_.is_solved && gui_ctrl_.is_init) { load_scene(); gui_ctrl_.is_init = false; } + viewer->data(0).show_lines = gui_ctrl_.mesh_visible; + load_bc(); + if (gui_ctrl_.is_modified) { redraw_scene(); gui_ctrl_.is_modified = false; @@ -196,7 +276,6 @@ namespace ipc { } - private: std::shared_ptr sp_StaticSim_; GUICtrl gui_ctrl_; diff --git a/static_sim/StaticSim.cpp b/static_sim/StaticSim.cpp index 31a0385..f221286 100644 --- a/static_sim/StaticSim.cpp +++ b/static_sim/StaticSim.cpp @@ -360,16 +360,18 @@ namespace ssim { } } - Model StaticSim::get_mesh() const { - Eigen::MatrixXd V_surf(SVI.size(), 3); - for (int svI = 0; svI < SVI.size(); ++svI) { - int vI = SVI(svI); - V_surf.row(svI) = TV.row(vI); + Model StaticSim::get_mesh() { + if (mesh_.NumVertex() == 0) { + // fill mesh_ + Eigen::MatrixXd V_surf(SVI.size(), 3); + for (int svI = 0; svI < SVI.size(); ++svI) { + int vI = SVI(svI); + V_surf.row(svI) = TV.row(vI); + } + mesh_.V = V_surf; + mesh_.F = F_surf; } - Model ret; - ret.V = V_surf; - ret.F = F_surf; - return ret; + return mesh_; } void StaticSim::MapAppendTarget(int target) { diff --git a/static_sim/StaticSim.h b/static_sim/StaticSim.h index 4d1f705..834e6db 100644 --- a/static_sim/StaticSim.h +++ b/static_sim/StaticSim.h @@ -14,14 +14,20 @@ #include "SimTargetOption.h" #ifdef LINSYSSOLVER_USE_CHOLMOD + #include "CHOLMODSolver.hpp" + #else #include "EigenLibSolver.hpp" #endif #include "BoundaryConditions.hpp" #include "Config.hpp" +#include +namespace ipc::rigid{ + class UIStaticSimState; +} namespace ssim { struct Model { Eigen::MatrixX3d V; @@ -36,112 +42,142 @@ namespace ssim { }; }; + class StaticSim { using MeshModel = Model; + using Mesh = Model; + friend class ipc::rigid::UIStaticSimState; public: - StaticSim(const SimTargetOption &option, const std::string &jsonPath); - - ~StaticSim() {} - - void simulation(); - - /** - * Given query points Q, compute their displacements and stress - * @param Q (nQ, 3), query points - * @param QU return value, (nQ), norm of displacements - * @param Qstress return value, (nQ, 6), stress size is (6) on each query point - */ - void postprocess(Eigen::MatrixXd &Q, - Eigen::VectorXd &QU, - Eigen::MatrixXd &Qstress); - - Eigen::MatrixXd EvaluateTarget(SimTargetOption::Target target) { - if (!option_.is_option_set(target)) { - // If no cache, update option_ and map_ - MapAppendTarget(target); - option_.set_option(target); - } - // Return cache - return map_target_to_evaluated_[target]; + + StaticSim(const SimTargetOption &option, const std::string &jsonPath); + + ~ + + StaticSim() {} + + void simulation(); + + + Mesh GetMeshInDBC() { + // TODO: get mesh in DBC box + Mesh ret; + ret.V = get_mesh().V; + ret.F = get_mesh().F.topRows(100); + + return ret; + } + + /** + * Given query points Q, compute their displacements and stress + * @param Q (nQ, 3), query points + * @param QU return value, (nQ), norm of displacements + * @param Qstress return value, (nQ, 6), stress size is (6) on each query point + */ + void postprocess(Eigen::MatrixXd &Q, + Eigen::VectorXd &QU, + Eigen::MatrixXd &Qstress); + + Eigen::MatrixXd EvaluateTarget(SimTargetOption::Target target) { + if (!option_.is_option_set(target)) { + // If no cache, update option_ and map_ + MapAppendTarget(target); + option_.set_option(target); } + // Return cache + return map_target_to_evaluated_[target]; + } - // return surf tri mesh of tet mesh - Model get_mesh() const; + // return surf tri mesh of tet mesh + Model get_mesh(); private: - void computeFeatures(); - void computeK(); - void solve(); - void setBC(); - void prepare_surf_result(); - - void MapAppendTarget(int target); - - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateUNorm() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateUX() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateUY() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateUZ() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSNorm() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSVonMises() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSX() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSY() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSZ() const; - // Return: 1 x 1 - Eigen::MatrixXd EvaluateCompliance() const; - private: - MeshModel mesh_; - SimTargetOption option_; - std::vector map_target_to_evaluated_; + void computeFeatures(); + + void computeK(); + + void solve(); + + void setBC(); + + void prepare_surf_result(); + + void MapAppendTarget(int target); + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUNorm() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUX() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUY() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUZ() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSNorm() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSVonMises() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSX() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSY() const; + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSZ() const; + + // Return: 1 x 1 + Eigen::MatrixXd EvaluateCompliance() const; private: - std::vector DirichletBCs; - std::vector NeumannBCs; - Eigen::Matrix D; // constitutive matrix - - // owned data - int nN, nEle, nDof; - Eigen::MatrixXd TV; // vertices coordinates - Eigen::MatrixXd TV1; // deformed vertices coordinates - Eigen::MatrixXi TT; // vertice index of each tetrahedron - Eigen::MatrixXi SF; - - int eleNodeNum; - int eleDofNum; - std::vector eDof; - - // owned features - Eigen::VectorXd load; // load of each dof - Eigen::VectorXd U; // dofs' displacement to be computed - - Eigen::VectorXi DBC_nI; // vertex in DBC - Eigen::VectorXi isDBC; // 0: not in DBC, 1: in DBC - - Eigen::VectorXi SVI; // vertice indices of surface nodes - Eigen::MatrixXi F_surf; // boundary vertice indices in surface triangles mesh - - // indices for fast access - std::vector> vNeighbor; // records all vertices' indices adjacent to each vertice - std::vector>> vFLoc; - std::shared_ptr> linSysSolver; - - // resulted data to evaluate - double compliance_; - Eigen::MatrixXd surf_U_; - Eigen::MatrixXd surf_stress_; - Eigen::VectorXd surf_vonstress_; + MeshModel mesh_; + SimTargetOption option_; + std::vector map_target_to_evaluated_; - }; + private: + std::vector DirichletBCs; + std::vector NeumannBCs; + Eigen::Matrix D; // constitutive matrix + + // owned data + int nN, nEle, nDof; + Eigen::MatrixXd TV; // vertices coordinates + Eigen::MatrixXd TV1; // deformed vertices coordinates + Eigen::MatrixXi TT; // vertice index of each tetrahedron + Eigen::MatrixXi SF; + + int eleNodeNum; + int eleDofNum; + std::vector eDof; + + // owned features + Eigen::VectorXd load; // load of each dof + Eigen::VectorXd U; // dofs' displacement to be computed + + Eigen::VectorXi DBC_nI; // vertex in DBC + Eigen::VectorXi isDBC; // 0: not in DBC, 1: in DBC + + Eigen::VectorXi SVI; // vertice indices of surface nodes + Eigen::MatrixXi F_surf; // boundary vertice indices in surface triangles mesh + + // indices for fast access + std::vector> vNeighbor; // records all vertices' indices adjacent to each vertice + std::vector>> vFLoc; + std::shared_ptr> linSysSolver; + + // resulted data to evaluate + double compliance_; + Eigen::MatrixXd surf_U_; + Eigen::MatrixXd surf_stress_; + Eigen::VectorXd surf_vonstress_; + +}; } // ssim