From 1864497497757a50694dfd47be4582605f97da80 Mon Sep 17 00:00:00 2001 From: Chenw831X Date: Fri, 14 Apr 2023 00:52:14 +0800 Subject: [PATCH] add color surface on BC --- src/viewer/UIStaticMenu.cpp | 56 +++++--- src/viewer/UIStaticSimState.h | 35 ++++- static_sim/StaticSim.cpp | 71 +++++++++- static_sim/StaticSim.h | 245 +++++++++++++++++++--------------- 4 files changed, 278 insertions(+), 129 deletions(-) diff --git a/src/viewer/UIStaticMenu.cpp b/src/viewer/UIStaticMenu.cpp index 11e3a37..6ab7519 100644 --- a/src/viewer/UIStaticMenu.cpp +++ b/src/viewer/UIStaticMenu.cpp @@ -90,7 +90,7 @@ namespace ipc::rigid { load(fname); } #else - load("/home/cflin/Documents/CppField/source/RigidElasticSim/sim-test/rigid-test/cube/config.json"); + load(CMAKE_SOURCE_DIR "/sim-test/rigid-test/rocker-arm/config.json"); #endif } if (m_has_scene) { @@ -150,14 +150,22 @@ namespace ipc::rigid { sp_StaticSim_->simulation(); gui_ctrl_.is_modified = true; gui_ctrl_.is_solved = true; + gui_ctrl_.single_color = false; + } + if (ImGui::Button(("原始模型"))) { + gui_ctrl_.single_color = true; + gui_ctrl_.is_modified = true; } // ImGui::RadioButton("开始计算##SimPlayer", &player_state, PlayerState::Playing); // ImGui::RadioButton("暂停##SimPlayer", &player_state, PlayerState::Paused); ImGui::Checkbox("网格可视化", &gui_ctrl_.mesh_visible); + // -------------------------------------------------------------------- if (ImGui::CollapsingHeader( "边界条件", ImGuiTreeNodeFlags_DefaultOpen)) { - ImGui::Checkbox("边界可视化", &gui_ctrl_.is_visible_BC); + if (ImGui::Checkbox("边界可视化", &gui_ctrl_.is_visible_BC)) { + gui_ctrl_.is_modified = true; + } if (ImGui::CollapsingHeader( "Dirichlet边界", ImGuiTreeNodeFlags_DefaultOpen)) { // -------------------------------------------------------------------- @@ -170,12 +178,19 @@ namespace ipc::rigid { (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::InputFloat3(("Dirichlet最小值 " + std::to_string(i + 1)).c_str(), v_min_point.data()); - ImGui::InputFloat3(("Dirichlet最大值 " + std::to_string(i + 1)).c_str(), v_max_point.data()); - 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]}; + if ( + ImGui::InputFloat3(("Dirichlet最小值 " + std::to_string(i + 1)).c_str(), + v_min_point.data()) || + ImGui::InputFloat3(("Dirichlet最大值 " + std::to_string(i + 1)).c_str(), + v_max_point.data())) { + gui_ctrl_.is_modified = true; + 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]}; + } } - sp_StaticSim_->updateBC(); + if (gui_ctrl_.is_modified) + sp_StaticSim_->updateBC(); + } @@ -190,16 +205,22 @@ namespace ipc::rigid { std::vector v_max_point = {(float) i_NBC.relMaxBBox.x(), (float) i_NBC.relMaxBBox.y(), (float) i_NBC.relMaxBBox.z()}; std::vector v_force = {(float) i_NBC.force.x(), (float) i_NBC.force.y(), - (float) i_NBC.force.z()}; - - ImGui::InputFloat3(("Neumann最小值 " + std::to_string(i + 1)).c_str(), v_min_point.data()); - ImGui::InputFloat3(("Neumann最大值 " + std::to_string(i + 1)).c_str(), v_max_point.data()); - ImGui::InputFloat3(("Neumann力 " + std::to_string(i + 1)).c_str(), v_force.data()); - i_NBC.relMinBBox = {v_min_point[0], v_min_point[1], v_min_point[2]}; - i_NBC.relMaxBBox = {v_max_point[0], v_max_point[1], v_max_point[2]}; - i_NBC.force={v_force[0],v_force[1],v_force[2]}; + (float) i_NBC.force.z()}; + if ( + ImGui::InputFloat3(("Neumann最小值 " + std::to_string(i + 1)).c_str(), + v_min_point.data()) || + ImGui::InputFloat3(("Neumann最大值 " + std::to_string(i + 1)).c_str(), + v_max_point.data()) || + ImGui::InputFloat3(("Neumann力 " + std::to_string(i + 1)).c_str(), v_force.data())) { + gui_ctrl_.is_modified = true; + i_NBC.relMinBBox = {v_min_point[0], v_min_point[1], v_min_point[2]}; + i_NBC.relMaxBBox = {v_max_point[0], v_max_point[1], v_max_point[2]}; + i_NBC.force = {v_force[0], v_force[1], v_force[2]}; + + } } - sp_StaticSim_->updateBC(); + if (gui_ctrl_.is_modified) + sp_StaticSim_->updateBC(); } @@ -254,7 +275,8 @@ namespace ipc::rigid { char *desc = "121.88"; // ImGui::Text("柔顺度(compliance)"); // sp_StaticSim_->EvaluateTarget(ssim::SimTargetOption::COMPLIANCE)(0,0); - ImGui::TextWrapped("柔顺度 (compliance): %.3e", sp_StaticSim_->EvaluateTarget(ssim::SimTargetOption::COMPLIANCE)(0,0)); + ImGui::TextWrapped("柔顺度 (compliance): %.3e", + sp_StaticSim_->EvaluateTarget(ssim::SimTargetOption::COMPLIANCE)(0, 0)); if (ImGui::IsItemHovered()) { // 显示工具的提示 ImGui::BeginTooltip(); diff --git a/src/viewer/UIStaticSimState.h b/src/viewer/UIStaticSimState.h index 6564d0c..5e30800 100644 --- a/src/viewer/UIStaticSimState.h +++ b/src/viewer/UIStaticSimState.h @@ -33,12 +33,13 @@ namespace ipc { bool is_loaded_json = false; bool is_visible_BC = false; bool mesh_visible = true; - bool single_color = false; + bool single_color = true; ssim::SimTargetOption::Target target_to_evaluate = ssim::SimTargetOption::U_NORM; void set_target(ssim::SimTargetOption::Target target_to_set) { target_to_evaluate = target_to_set; is_modified = true; + single_color = false; } }; @@ -141,6 +142,8 @@ namespace ipc { ssim::Model model = sp_StaticSim_->get_mesh(); const Eigen::MatrixXd &V = model.V; const Eigen::MatrixXi &F = model.F; + Eigen::VectorXi DBC_vertex_idx = sp_StaticSim_->DBC_vertexIdx_; + Eigen::VectorXi NBC_vertex_idx = sp_StaticSim_->NBC_vertexIdx_; Eigen::MatrixXd target = sp_StaticSim_->EvaluateTarget(gui_ctrl_.target_to_evaluate); @@ -155,10 +158,30 @@ namespace ipc { // colormap Eigen::MatrixXd C; igl::jet(target, true, C); + if (0 && gui_ctrl_.is_visible_BC) { + for (int i_ = 0; i_ < DBC_vertex_idx.size(); ++i_) { + C.row(DBC_vertex_idx(i_)) = Eigen::RowVector3d(0, 0, 1); + } + for (int i_ = 0; i_ < NBC_vertex_idx.size(); ++i_) { + C.row(NBC_vertex_idx(i_)) = Eigen::RowVector3d(1, 0, 0); + } + } viewer->data(0).set_colors(C); } else { + Eigen::MatrixXd color(V.rows(), 3); + for (int i_ = 0; i_ < V.rows(); ++i_) { + color.row(i_) = Eigen::RowVector3d(1, 1, 0); + } + if (gui_ctrl_.is_visible_BC) { + for (int i_ = 0; i_ < DBC_vertex_idx.size(); ++i_) { + color.row(DBC_vertex_idx(i_)) = Eigen::RowVector3d(0, 0, 1); + } + for (int i_ = 0; i_ < NBC_vertex_idx.size(); ++i_) { + color.row(NBC_vertex_idx(i_)) = Eigen::RowVector3d(1, 0, 0); + } + } // single color - viewer->data(0).set_colors(Eigen::RowVector3d(1, 0, 0)); + viewer->data(0).set_colors(color); // colorbar_plugin_.draw_colorbar(igl::ColorMapType::COLOR_MAP_TYPE_JET,0,1,/*background_color=*/Eigen::Vector4f(1,1,1,0)); } @@ -220,17 +243,21 @@ namespace ipc { for (int i = 0; i < sp_StaticSim_->DirichletBCs.size(); ++i) { auto [V, F] = get_V_F_from_bd_box(sp_StaticSim_->DirichletBCs[i].absMinBBox, sp_StaticSim_->DirichletBCs[i].absMaxBBox); + viewer->data(i + 1).clear(); + viewer->data(i + 1).show_lines = false; viewer->data(i + 1).set_mesh(V, F); - viewer->data(i + 1).set_colors(Eigen::RowVector4d(1, 0, 0, 0.5)); + viewer->data(i + 1).set_colors(Eigen::RowVector4d(0, 0, 0, 0.1)); 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].absMinBBox, sp_StaticSim_->NeumannBCs[i].absMaxBBox); + viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).clear(); + viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).show_lines = false; viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).set_mesh(V, F); viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).set_colors( - Eigen::RowVector4d(0, 0, 1, 0.5)); + Eigen::RowVector4d(0, 0, 0, 0.1)); viewer->data(i + 1 + sp_StaticSim_->DirichletBCs.size()).is_visible = gui_ctrl_.is_visible_BC; } diff --git a/static_sim/StaticSim.cpp b/static_sim/StaticSim.cpp index 028ef9b..5562175 100644 --- a/static_sim/StaticSim.cpp +++ b/static_sim/StaticSim.cpp @@ -25,6 +25,9 @@ namespace ssim { exit(-1); } + int ENUM_SIZE = SimTargetOption::Target::ENUM_SIZE; + map_target_to_evaluated_.resize(ENUM_SIZE); + DirichletBCs = config.DirichletBCs; NeumannBCs = config.NeumannBCs; @@ -106,6 +109,73 @@ namespace ssim { for (auto &NBC: NeumannBCs) { NBC.calcAbsBBox(modelMinBBox, modelMaxBBox); } + + { + int cnt = 0; + for (auto &DBC: DirichletBCs) { + Eigen::VectorXi DBC_i_faceIdx = getSurfTriForBox(DBC.absMinBBox, DBC.absMaxBBox); + DBC_faceIdx_.conservativeResize(cnt + DBC_i_faceIdx.size()); + DBC_faceIdx_.segment(cnt, DBC_i_faceIdx.size()) = DBC_i_faceIdx; + } + std::set DBC_faceIdx_set(DBC_faceIdx_.data(), DBC_faceIdx_.data() + DBC_faceIdx_.size()); + std::vector DBC_faceIdx_vec(DBC_faceIdx_set.begin(), DBC_faceIdx_set.end()); + DBC_faceIdx_ = Eigen::Map(DBC_faceIdx_vec.data(), DBC_faceIdx_vec.size());; + + std::set DBC_vertex_set; + for (int fI_ = 0; fI_ < DBC_faceIdx_.size(); ++fI_) { + int fI = DBC_faceIdx_(fI_); + DBC_vertex_set.insert(vI2SVI[SF(fI, 0)]); + DBC_vertex_set.insert(vI2SVI[SF(fI, 1)]); + DBC_vertex_set.insert(vI2SVI[SF(fI, 2)]); + } + std::vector DBC_vertexIdx_vec(DBC_vertex_set.begin(), DBC_vertex_set.end()); + DBC_vertexIdx_ = Eigen::Map(DBC_vertexIdx_vec.data(), DBC_vertexIdx_vec.size());; + } + { + int cnt = 0; + for (auto &NBC: NeumannBCs) { + Eigen::VectorXi NBC_i_faceIdx = getSurfTriForBox(NBC.absMinBBox, NBC.absMaxBBox); + NBC_faceIdx_.conservativeResize(cnt + NBC_i_faceIdx.size()); + NBC_faceIdx_.segment(cnt, NBC_i_faceIdx.size()) = NBC_i_faceIdx; + } + std::set NBC_faceIdx_set(NBC_faceIdx_.data(), NBC_faceIdx_.data() + NBC_faceIdx_.size()); + std::vector NBC_faceIdx_vec(NBC_faceIdx_set.begin(), NBC_faceIdx_set.end()); + NBC_faceIdx_ = Eigen::Map(NBC_faceIdx_vec.data(), NBC_faceIdx_vec.size());; + + std::set NBC_vertex_set; + for (int fI_ = 0; fI_ < NBC_faceIdx_.size(); ++fI_) { + int fI = NBC_faceIdx_(fI_); + NBC_vertex_set.insert(vI2SVI[SF(fI, 0)]); + NBC_vertex_set.insert(vI2SVI[SF(fI, 1)]); + NBC_vertex_set.insert(vI2SVI[SF(fI, 2)]); + } + std::vector NBC_vertexIdx_vec(NBC_vertex_set.begin(), NBC_vertex_set.end()); + NBC_vertexIdx_ = Eigen::Map(NBC_vertexIdx_vec.data(), NBC_vertexIdx_vec.size());; + } + } + + Eigen::VectorXi StaticSim::getSurfTriForBox(const Eigen::Vector3d &minBox, + const Eigen::Vector3d &maxBox) { + Eigen::VectorXi face_Idx(SF.rows()); + + int f_cnt = 0; + for (int sfI = 0; sfI < SF.rows(); ++sfI) { + bool flag = true; + for (int i_ = 0; i_ < 3; ++i_) { + const Eigen::Vector3d &P = TV.row(SF(sfI, i_)); + if (((P - minBox).array() < 0).any() || ((P - maxBox).array() > 0).any()) { + flag = false; + break; + } + } + + if (flag) { + face_Idx(f_cnt++) = sfI; + } + } + + face_Idx.conservativeResize(f_cnt); + return face_Idx; } void StaticSim::prepare_surf_result() { @@ -172,7 +242,6 @@ namespace ssim { void StaticSim::computeFeatures() { // compute F_surf int cnt = 0; - std::unordered_map vI2SVI; for (int sfI = 0; sfI < SF.rows(); ++sfI) { for (int j = 0; j < 3; ++j) { const int &vI = SF(sfI, j); diff --git a/static_sim/StaticSim.h b/static_sim/StaticSim.h index b85b598..cc15ab4 100644 --- a/static_sim/StaticSim.h +++ b/static_sim/StaticSim.h @@ -8,6 +8,8 @@ #include #include #include +#include +#include #include #include #include @@ -25,7 +27,7 @@ #include "Config.hpp" #include -namespace ipc::rigid{ +namespace ipc::rigid { class UIStaticSimState; } namespace ssim { @@ -42,149 +44,178 @@ namespace ssim { }; }; - struct MaterialProperty{ - float Youngs_Modulus=1.0f; - float Poisson_ratio=0.3f; - float density=1.0f; + struct MaterialProperty { + float Youngs_Modulus = 1.0f; + float Poisson_ratio = 0.3f; + float density = 1.0f; }; 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(); - - /** - * update BC's absBBox - */ - void updateBC(); - - - /** - * 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); + StaticSim(const SimTargetOption &option, const std::string &jsonPath); + + ~StaticSim() {} + + void simulation(); + + /** + * update BC's absBBox + */ + void updateBC(); + + /** + * 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) || map_target_to_evaluated_[target].size()==0) { + // If no cache, update option_ and map_ + MapAppendTarget(target); + option_.set_option(target); + } + // Return cache + return map_target_to_evaluated_[target]; } - // Return cache - return map_target_to_evaluated_[target]; - } - // return surf tri mesh of tet mesh - Model get_mesh(); + // return surf tri mesh of tet mesh + Model get_mesh(); private: - void computeFeatures(); + void computeFeatures(); - void computeK(); + void computeK(); - void solve(); + void solve(); - void setBC(); + void setBC(); - void prepare_surf_result(); + void prepare_surf_result(); - void MapAppendTarget(int target); + void MapAppendTarget(int target); - MaterialProperty& get_material_property(){ - return material_property_; - } + Eigen::VectorXi getSurfTriForBox(const Eigen::Vector3d &minBox, + const Eigen::Vector3d &maxBox); + MaterialProperty &get_material_property() { + return material_property_; + } + + + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUNorm() const; - // 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 EvaluateUX() const; + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateUY() 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 EvaluateUZ() const; + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSNorm() 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 EvaluateSVonMises() const; + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSX() 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 EvaluateSY() const; + // Return: #mesh.V.rows() x 1 + Eigen::MatrixXd EvaluateSZ() const; - // Return: #mesh.V.rows() x 1 - Eigen::MatrixXd EvaluateSZ() const; + // Return: 1 x 1 + Eigen::MatrixXd EvaluateCompliance() const; - // Return: 1 x 1 - Eigen::MatrixXd EvaluateCompliance() const; + void writePntVTK(const std::string &path, const Eigen::MatrixXd &V) { + std::ofstream out(path); + out << "# vtk DataFile Version 3.0\n" + "Volume Mesh\n" + "ASCII\n" + "DATASET UNSTRUCTURED_GRID" << std::endl; + out << "POINTS " << V.rows() << " float" << std::endl; + for (int i = 0; i < V.rows(); ++i) { + out << std::setprecision(4) << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << std::endl; + } + out << "CELLS " << V.rows() << " " << V.rows() * (1 + 1) << std::endl; + for (int i = 0; i < V.rows(); ++i) { + out << "1 " << i << std::endl; + } + out << "CELL_TYPES " << V.rows() << std::endl; + for (int i = 0; i < V.rows(); ++i) { + out << 1 << std::endl; + } + } private: - MeshModel mesh_; - MaterialProperty material_property_; - SimTargetOption option_; - std::vector map_target_to_evaluated_; + MeshModel mesh_; + MaterialProperty material_property_; + 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_; - -}; + 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 + std::unordered_map vI2SVI; + + // 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_; + + Eigen::VectorXi DBC_faceIdx_; + Eigen::VectorXi DBC_vertexIdx_; + Eigen::VectorXi NBC_faceIdx_; + Eigen::VectorXi NBC_vertexIdx_; + + }; } // ssim