diff --git a/README.md b/README.md index cee79a2..660f64c 100644 --- a/README.md +++ b/README.md @@ -70,6 +70,16 @@ make -j 16 ``` ## Usage +### 24/5/8 updated: extract sequence vtk for topo-optimization +Usage: +1. set `extract_inner_vtk` to `true` in `example/top_mechanical/main.cpp` or `example/top_thermoelastic/main.cpp`. +```c++ +// example/top_mechanical/main.cpp or +// example/top_thermoelastic/main.cpp + bool extract_inner_vtk = true; +``` +2. run `example/top_mechanical` or `example/top_thermoelastic`. +3. see `docs/paraview_extract_animation.webm` for how to extract animation in Paraview. ### updated: top and sim sequence 1. Use `example/defined_model_writer` to generate a user-defined voxel mesh/model as TO input. 2. Use `example/top_mechanical` to run mechanical TO. diff --git a/docs/paraview_extract_animation.webm b/docs/paraview_extract_animation.webm new file mode 100644 index 0000000..a8a66ef Binary files /dev/null and b/docs/paraview_extract_animation.webm differ diff --git a/examples/top_mechanical/main.cpp b/examples/top_mechanical/main.cpp index e9e54e9..a4175c8 100644 --- a/examples/top_mechanical/main.cpp +++ b/examples/top_mechanical/main.cpp @@ -6,6 +6,7 @@ int main() { fs_path config_file(ASSETS_DIR "/top/config_Lshape.json"); fs_path model_file( ASSETS_DIR "/voxel_model/Lshape_model.txt"); + bool extract_inner_vtk = true; // END INPUT std::string model_name = model_file.filename().replace_extension(); @@ -24,7 +25,7 @@ int main() { auto sp_mech_top3d = sp_mech_inter->CreatTop(); // loop spdlog::critical("start to mechanical top opt ..."); - Tensor3d t_me_rho = sp_mech_top3d->TopOptMainLoop(); + Tensor3d t_me_rho = sp_mech_top3d->TopOptMainLoop(false, extract_inner_vtk); // postprocess { spdlog::critical("extract compliance and volume each iteration ..."); @@ -54,4 +55,23 @@ int main() { WriteTensorToVtk(rho_vtk_path, t_me_rho, sp_mech_inter->sp_mesh_); spdlog::info("write density vtk to: {}", rho_vtk_path.c_str()); } + if (extract_inner_vtk) { + if (sp_mech_top3d->v_ten_rho_.size() == 0) { + throw std::runtime_error("v_ten_rho_ is empty"); + } + for (int i = 0; i v_ten_rho_.size(); ++i) { + std::stringstream ss; + ss << std::setw(3) << std::setfill('0') << i; + std::string suffix_idx = ss.str(); + fs_path inner_vtk_path = + output_dir / "vtk" / ex_name / + (ex_name + "_MeTop" + "_inner") / (suffix_idx + ".vtk"); + WriteTensorToVtk(inner_vtk_path,sp_mech_top3d->v_ten_rho_[i], + sp_mech_inter->sp_mesh_); + if (i == 0) { + spdlog::info("write inner vtk(0,1,...) to: {}", + inner_vtk_path.c_str()); + } + } + } } \ No newline at end of file diff --git a/examples/top_thermoelastic/main.cpp b/examples/top_thermoelastic/main.cpp index 95ae9cd..ae94333 100644 --- a/examples/top_thermoelastic/main.cpp +++ b/examples/top_thermoelastic/main.cpp @@ -1,11 +1,16 @@ +#include #include "ConfigMechanicalInterface.h" #include "ConfigThermalInterface.h" int main() { // INPUT ARGS - fs_path config_file(ASSETS_DIR "/top/config_Lshape.json"); + fs_path + config_file(ASSETS_DIR + "/top/config_Lshape.json"); fs_path model_file( - ASSETS_DIR "/voxel_model/Lshape_model.txt"); + ASSETS_DIR + "/voxel_model/Lshape_model.txt"); + bool extract_inner_vtk = true; // END INPUT std::string model_name = model_file.filename().replace_extension(); @@ -15,10 +20,10 @@ int main() { spdlog::info("Read model from '{}'", model_file.string()); spdlog::info("Output to '{}'", output_dir.string()); - std::shared_ptr sp_mech_inter = std::make_shared( - config_file,tr_model); - std::shared_ptr sp_ther_inter = std::make_shared( - config_file,tr_model); + std::shared_ptr sp_mech_inter = std::make_shared( + config_file, tr_model); + std::shared_ptr sp_ther_inter = std::make_shared( + config_file, tr_model); std::string ex_name = sp_mech_inter->ex_name_; spdlog::critical("Thermoelastic TO example: {}", ex_name); @@ -30,7 +35,8 @@ int main() { // init thermoelastic top3d top::ThermoelasticTop3d mech_ther_top3d(sp_mech_top3d, sp_ther_top3d); // loop - top::Tensor3d t_meth_rho = mech_ther_top3d.TopOptMainLoop(); + top::Tensor3d t_meth_rho = mech_ther_top3d.TopOptMainLoop(false, + extract_inner_vtk); // postprocess { spdlog::critical("extract compliance and volume each iteration ..."); @@ -60,4 +66,18 @@ int main() { WriteTensorToVtk(rho_vtk_path, t_meth_rho, sp_mech_inter->sp_mesh_); spdlog::info("write density vtk to: {}", rho_vtk_path.c_str()); } + if (extract_inner_vtk) { + for (int i = 0; i < sp_ther_top3d->v_ten_rho_.size(); ++i) { + std::stringstream ss; + ss << std::setw(3) << std::setfill('0') << i; + fs_path inner_vtk_path = + output_dir / "vtk" / ex_name / + (ex_name + "_MeThTop" + "_inner") / ss.str() / ".vtk"; + WriteTensorToVtk(inner_vtk_path, t_meth_rho, + sp_mech_inter->sp_mesh_, true); + if (i == 0) { + spdlog::info("write inner vtk(0,1,...) to: {}", inner_vtk_path.c_str()); + } + } + } } \ No newline at end of file diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index f061ea1..0cd6620 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -1,3 +1,4 @@ + // // Created by cflin on 6/12/23. // @@ -6,15 +7,18 @@ #include "LinearSolver/Amgcl.h" #ifdef USE_AMGCL_CUDA + #include "LinearSolver/AmgclCuda.h" + #endif da::sha::top::Tensor3d -da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { +da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation, + bool save_internal_rho) { auto &sp_mesh_ = sp_mech_top3d_->sp_mesh_; auto &sp_para_ = sp_mech_top3d_->sp_para_; Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); - xPhys_col.setConstant(only_simulation?1:sp_para_->volfrac); + xPhys_col.setConstant(only_simulation ? 1 : sp_para_->volfrac); bool flg_chosen = false; Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); // bool flg_chosen = chosen_ele_id.size() != 0; @@ -42,9 +46,6 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { dv.setOnes(); dv = (sp_mech_top3d_->H_ * dv).array() / sp_mech_top3d_->Hs_.array(); - Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); Eigen::VectorXi global_idx_of_ele_in_use = sp_mesh_->GetGlobalIdxOfEleInUse(); // dofs of limited T @@ -83,6 +84,36 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { LOG_SOLVER + // set 0 to rho of unchosen part + assert(xPhys_col.size()); + Eigen::VectorXi continue_idx = + Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, + chosen_ele_id) + : Eigen::VectorXi(); + + auto RhoVec2Ten = [&](Eigen::VectorXd rho) -> Tensor3d { + rho(unchosen_idx).setZero(); + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz()); + ele_to_write(global_idx_of_ele_in_use) = rho; + Tensor3d ten_xPhys_to_write( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, + 1); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape( + Eigen::array{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); + return ten_xPhys_to_write; + }; + + if (save_internal_rho) + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); // start iteration while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; @@ -348,9 +379,9 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { #define SENSITIVITY_SCALE_COEF 200 size_t num_constraints = 1 - #ifndef MECH_ONLY - #ifdef WITH_T_LIMIT - + dT_dx.cols()// volume and temperature constraints +#ifndef MECH_ONLY +#ifdef WITH_T_LIMIT + + dT_dx.cols()// volume and temperature constraints #endif #endif ; @@ -384,7 +415,7 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { << 1.0 / num_variables * dv #ifndef MECH_ONLY #ifdef WITH_T_LIMIT - , 1.0 / sp_para_->T_limit * dT_dx + , 1.0 / sp_para_->T_limit * dT_dx #endif #endif ).finished().transpose() * @@ -409,55 +440,20 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { } spdlog::info("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", - loop, c, v, - change); + loop, c, v, + change); v_compliance_.push_back(c); v_volume_.push_back(v); if (only_simulation) { break; } + if(save_internal_rho) + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); } // result sp_mech_top3d_->rho_ = only_simulation ? sp_mesh_->GetInitEleRho() : xPhys_col; - // set 0 to rho of unchosen part - assert(xPhys_col.size()); - Eigen::VectorXi continue_idx = - Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, - chosen_ele_id) - : Eigen::VectorXi(); - { - xPhys_col(unchosen_idx).setZero(); - ele_to_write(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - sp_mech_top3d_->rho_field_zero_filled_ = ten_xPhys_to_write; - } - - { - xPhys_col(unchosen_idx).setOnes(); - ele_to_write(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - sp_mech_top3d_->rho_field_one_filled_ = ten_xPhys_to_write; - } + sp_mech_top3d_->rho_field_zero_filled_ = RhoVec2Ten(xPhys_col); return sp_mech_top3d_->rho_field_zero_filled_; } diff --git a/src/ThermoelasticTop3d.cu b/src/ThermoelasticTop3d.cu index f061ea1..0cd6620 100644 --- a/src/ThermoelasticTop3d.cu +++ b/src/ThermoelasticTop3d.cu @@ -1,3 +1,4 @@ + // // Created by cflin on 6/12/23. // @@ -6,15 +7,18 @@ #include "LinearSolver/Amgcl.h" #ifdef USE_AMGCL_CUDA + #include "LinearSolver/AmgclCuda.h" + #endif da::sha::top::Tensor3d -da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { +da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation, + bool save_internal_rho) { auto &sp_mesh_ = sp_mech_top3d_->sp_mesh_; auto &sp_para_ = sp_mech_top3d_->sp_para_; Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); - xPhys_col.setConstant(only_simulation?1:sp_para_->volfrac); + xPhys_col.setConstant(only_simulation ? 1 : sp_para_->volfrac); bool flg_chosen = false; Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); // bool flg_chosen = chosen_ele_id.size() != 0; @@ -42,9 +46,6 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { dv.setOnes(); dv = (sp_mech_top3d_->H_ * dv).array() / sp_mech_top3d_->Hs_.array(); - Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); Eigen::VectorXi global_idx_of_ele_in_use = sp_mesh_->GetGlobalIdxOfEleInUse(); // dofs of limited T @@ -83,6 +84,36 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { LOG_SOLVER + // set 0 to rho of unchosen part + assert(xPhys_col.size()); + Eigen::VectorXi continue_idx = + Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, + chosen_ele_id) + : Eigen::VectorXi(); + + auto RhoVec2Ten = [&](Eigen::VectorXd rho) -> Tensor3d { + rho(unchosen_idx).setZero(); + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz()); + ele_to_write(global_idx_of_ele_in_use) = rho; + Tensor3d ten_xPhys_to_write( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, + 1); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape( + Eigen::array{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); + return ten_xPhys_to_write; + }; + + if (save_internal_rho) + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); // start iteration while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; @@ -348,9 +379,9 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { #define SENSITIVITY_SCALE_COEF 200 size_t num_constraints = 1 - #ifndef MECH_ONLY - #ifdef WITH_T_LIMIT - + dT_dx.cols()// volume and temperature constraints +#ifndef MECH_ONLY +#ifdef WITH_T_LIMIT + + dT_dx.cols()// volume and temperature constraints #endif #endif ; @@ -384,7 +415,7 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { << 1.0 / num_variables * dv #ifndef MECH_ONLY #ifdef WITH_T_LIMIT - , 1.0 / sp_para_->T_limit * dT_dx + , 1.0 / sp_para_->T_limit * dT_dx #endif #endif ).finished().transpose() * @@ -409,55 +440,20 @@ da::sha::top::ThermoelasticTop3d::TopOptMainLoop(bool only_simulation) { } spdlog::info("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", - loop, c, v, - change); + loop, c, v, + change); v_compliance_.push_back(c); v_volume_.push_back(v); if (only_simulation) { break; } + if(save_internal_rho) + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); } // result sp_mech_top3d_->rho_ = only_simulation ? sp_mesh_->GetInitEleRho() : xPhys_col; - // set 0 to rho of unchosen part - assert(xPhys_col.size()); - Eigen::VectorXi continue_idx = - Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, - chosen_ele_id) - : Eigen::VectorXi(); - { - xPhys_col(unchosen_idx).setZero(); - ele_to_write(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - sp_mech_top3d_->rho_field_zero_filled_ = ten_xPhys_to_write; - } - - { - xPhys_col(unchosen_idx).setOnes(); - ele_to_write(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - sp_mech_top3d_->rho_field_one_filled_ = ten_xPhys_to_write; - } + sp_mech_top3d_->rho_field_zero_filled_ = RhoVec2Ten(xPhys_col); return sp_mech_top3d_->rho_field_zero_filled_; } diff --git a/src/ThermoelasticTop3d.h b/src/ThermoelasticTop3d.h index 8d10445..3ca3188 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -17,13 +17,14 @@ namespace da::sha::top { Precompute(); } - Tensor3d TopOptMainLoop(bool only_simulation=false); + Tensor3d TopOptMainLoop(bool only_simulation = false, + bool save_internal_rho = false); Eigen::VectorXd GetTemperature() const { return sp_thermal_top3d_->GetU(); } - Eigen::VectorXd GetNormedDisplacement()const{ + Eigen::VectorXd GetNormedDisplacement() const { return sp_mech_top3d_->GetNormedDisplacement(); }; @@ -70,6 +71,7 @@ namespace da::sha::top { public: std::vector v_compliance_; std::vector v_volume_; + std::vector v_ten_rho_; private: std::shared_ptr sp_mech_top3d_, sp_thermal_top3d_; diff --git a/src/Top3d.cpp b/src/Top3d.cpp index 42ec804..4edd8b7 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -14,7 +14,8 @@ #endif namespace da::sha { namespace top { - Tensor3d Top3d::TopOptMainLoop(bool only_simulation) { + Tensor3d + Top3d::TopOptMainLoop(bool only_simulation, bool save_internal_rho) { Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); bool flg_chosen = chosen_ele_id.size() != 0; @@ -49,6 +50,40 @@ namespace da::sha { LOG_SOLVER + // set 0 to rho of unchosen part + assert(xPhys_col.size()); + Eigen::VectorXi continue_idx = + Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, + xPhys_col.size() - 1); + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference( + continue_idx, chosen_ele_id) + : Eigen::VectorXi(); + + auto RhoVec2Ten = [&](Eigen::VectorXd rho) -> Tensor3d { + rho(unchosen_idx).setZero(); + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz()); + ele_to_write(pixel_idx) = rho; + Tensor3d ten_xPhys_to_write( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz(), 1, + 1); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape( + Eigen::array < Eigen::DenseIndex, 3 > { + sp_mesh_->GetLx(), sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); + return ten_xPhys_to_write; + }; + + if (save_internal_rho) { + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); + } + // start iteration while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; @@ -128,54 +163,18 @@ namespace da::sha { if (only_simulation) { break; } + if (save_internal_rho) { + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); + } } // result rho_ = only_simulation ? sp_mesh_->GetInitEleRho() : xPhys_col; - // set 0 to rho of unchosen part - assert(xPhys_col.size()); - Eigen::VectorXi continue_idx = - Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, - xPhys_col.size() - 1); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference( - continue_idx, chosen_ele_id) - : Eigen::VectorXi(); - { - xPhys_col(unchosen_idx).setZero(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * - sp_mesh_->GetLz(), 1, 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - rho_field_zero_filled_ = ten_xPhys_to_write; - } - - { - xPhys_col(unchosen_idx).setOnes(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * - sp_mesh_->GetLz(), 1, 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - rho_field_one_filled_ = ten_xPhys_to_write; - } - return rho_field_zero_filled_; + return RhoVec2Ten(rho_); } - std::vector + std::vector Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) { Eigen::VectorXd ele_to_write = Eigen::VectorXd::Zero( @@ -191,7 +190,7 @@ namespace da::sha { mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue; } // fill - std::vector vt; + std::vector vt; for (int i = 0; i < which_col_of_stress.size(); ++i) { ele_to_write(pixel_idx) = mat_stress.col( which_col_of_stress(i)); @@ -208,9 +207,10 @@ namespace da::sha { sp_mesh_->GetLz()); Eigen::VectorXi pixel_idx = sp_mesh_->GetGlobalIdxOfEleInUse(); // stress - Eigen::VectorXd node_to_write=Eigen::VectorXd::Zero(sp_mesh_->GetNumNodes()); + Eigen::VectorXd node_to_write = Eigen::VectorXd::Zero( + sp_mesh_->GetNumNodes()); - std::vector v_B(8); + std::vector v_B(8); static const Eigen::MatrixXd delta_coord = (Eigen::MatrixXd(8, 3) << 0, 0, 0, 1, 0, 0, @@ -221,8 +221,8 @@ namespace da::sha { 1, 1, 1, 0, 1, 1 ).finished(); - for(int i=0;icomputeBe(delta_coord.row(i).array()-0.5); + for (int i = 0; i < delta_coord.rows(); ++i) { + v_B[i] = sp_fea_->computeBe(delta_coord.row(i).array() - 0.5); } auto computeVonStress = [&](Eigen::VectorXd stress) -> double { double x = stress(0); @@ -235,17 +235,19 @@ namespace da::sha { 3 * (xy * xy + yz * yz + zx * zx)); }; for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i);// 1x 24 + Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs( + i);// 1x 24 Eigen::VectorXd Ue = U_(dofs_in_ele_i); - Eigen::VectorXi node_id_in_ele_i = sp_mesh_->MapEleId2NodeIds(i);// 1x8 - for(int j=0;jMapEleId2NodeIds( + i);// 1x8 + for (int j = 0; j < node_id_in_ele_i.size(); ++j) { + int i_node = node_id_in_ele_i(j); + if (node_to_write(i_node) != 0) { continue; } - Eigen::MatrixXd B=v_B[j]; - Eigen::VectorXd s6=rho_(i) * sp_fea_->computeD() * B * Ue; - node_to_write(i_node)=computeVonStress(s6); + Eigen::MatrixXd B = v_B[j]; + Eigen::VectorXd s6 = rho_(i) * sp_fea_->computeD() * B * Ue; + node_to_write(i_node) = computeVonStress(s6); } } @@ -262,9 +264,9 @@ namespace da::sha { ten_prop_to_write(i, 0, 0) = proprty_col(i); } ten_prop_to_write = ten_prop_to_write.reshape( - Eigen::array{sp_mesh_->GetLx(), - sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); + Eigen::array < Eigen::DenseIndex, 3 > {sp_mesh_->GetLx(), + sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); return ten_prop_to_write; } @@ -339,7 +341,8 @@ namespace da::sha { } } } - std::vector> v_tri = Vec2Triplet(iH, jH, sH); + std::vector > v_tri = Vec2Triplet(iH, jH, + sH); H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles()); H_.setFromTriplets(v_tri.begin(), v_tri.end()); diff --git a/src/Top3d.cu b/src/Top3d.cu index 42ec804..4edd8b7 100644 --- a/src/Top3d.cu +++ b/src/Top3d.cu @@ -14,7 +14,8 @@ #endif namespace da::sha { namespace top { - Tensor3d Top3d::TopOptMainLoop(bool only_simulation) { + Tensor3d + Top3d::TopOptMainLoop(bool only_simulation, bool save_internal_rho) { Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); bool flg_chosen = chosen_ele_id.size() != 0; @@ -49,6 +50,40 @@ namespace da::sha { LOG_SOLVER + // set 0 to rho of unchosen part + assert(xPhys_col.size()); + Eigen::VectorXi continue_idx = + Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, + xPhys_col.size() - 1); + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference( + continue_idx, chosen_ele_id) + : Eigen::VectorXi(); + + auto RhoVec2Ten = [&](Eigen::VectorXd rho) -> Tensor3d { + rho(unchosen_idx).setZero(); + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz()); + ele_to_write(pixel_idx) = rho; + Tensor3d ten_xPhys_to_write( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * + sp_mesh_->GetLz(), 1, + 1); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape( + Eigen::array < Eigen::DenseIndex, 3 > { + sp_mesh_->GetLx(), sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); + return ten_xPhys_to_write; + }; + + if (save_internal_rho) { + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); + } + // start iteration while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; @@ -128,54 +163,18 @@ namespace da::sha { if (only_simulation) { break; } + if (save_internal_rho) { + v_ten_rho_.emplace_back(RhoVec2Ten(xPhys_col)); + } } // result rho_ = only_simulation ? sp_mesh_->GetInitEleRho() : xPhys_col; - // set 0 to rho of unchosen part - assert(xPhys_col.size()); - Eigen::VectorXi continue_idx = - Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, - xPhys_col.size() - 1); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference( - continue_idx, chosen_ele_id) - : Eigen::VectorXi(); - { - xPhys_col(unchosen_idx).setZero(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * - sp_mesh_->GetLz(), 1, 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - rho_field_zero_filled_ = ten_xPhys_to_write; - } - - { - xPhys_col(unchosen_idx).setOnes(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write( - sp_mesh_->GetLx() * sp_mesh_->GetLy() * - sp_mesh_->GetLz(), 1, 1); - for (int i = 0; i < ele_to_write.size(); ++i) { - ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); - } - ten_xPhys_to_write = ten_xPhys_to_write.reshape( - Eigen::array{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); - rho_field_one_filled_ = ten_xPhys_to_write; - } - return rho_field_zero_filled_; + return RhoVec2Ten(rho_); } - std::vector + std::vector Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) { Eigen::VectorXd ele_to_write = Eigen::VectorXd::Zero( @@ -191,7 +190,7 @@ namespace da::sha { mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue; } // fill - std::vector vt; + std::vector vt; for (int i = 0; i < which_col_of_stress.size(); ++i) { ele_to_write(pixel_idx) = mat_stress.col( which_col_of_stress(i)); @@ -208,9 +207,10 @@ namespace da::sha { sp_mesh_->GetLz()); Eigen::VectorXi pixel_idx = sp_mesh_->GetGlobalIdxOfEleInUse(); // stress - Eigen::VectorXd node_to_write=Eigen::VectorXd::Zero(sp_mesh_->GetNumNodes()); + Eigen::VectorXd node_to_write = Eigen::VectorXd::Zero( + sp_mesh_->GetNumNodes()); - std::vector v_B(8); + std::vector v_B(8); static const Eigen::MatrixXd delta_coord = (Eigen::MatrixXd(8, 3) << 0, 0, 0, 1, 0, 0, @@ -221,8 +221,8 @@ namespace da::sha { 1, 1, 1, 0, 1, 1 ).finished(); - for(int i=0;icomputeBe(delta_coord.row(i).array()-0.5); + for (int i = 0; i < delta_coord.rows(); ++i) { + v_B[i] = sp_fea_->computeBe(delta_coord.row(i).array() - 0.5); } auto computeVonStress = [&](Eigen::VectorXd stress) -> double { double x = stress(0); @@ -235,17 +235,19 @@ namespace da::sha { 3 * (xy * xy + yz * yz + zx * zx)); }; for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i);// 1x 24 + Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs( + i);// 1x 24 Eigen::VectorXd Ue = U_(dofs_in_ele_i); - Eigen::VectorXi node_id_in_ele_i = sp_mesh_->MapEleId2NodeIds(i);// 1x8 - for(int j=0;jMapEleId2NodeIds( + i);// 1x8 + for (int j = 0; j < node_id_in_ele_i.size(); ++j) { + int i_node = node_id_in_ele_i(j); + if (node_to_write(i_node) != 0) { continue; } - Eigen::MatrixXd B=v_B[j]; - Eigen::VectorXd s6=rho_(i) * sp_fea_->computeD() * B * Ue; - node_to_write(i_node)=computeVonStress(s6); + Eigen::MatrixXd B = v_B[j]; + Eigen::VectorXd s6 = rho_(i) * sp_fea_->computeD() * B * Ue; + node_to_write(i_node) = computeVonStress(s6); } } @@ -262,9 +264,9 @@ namespace da::sha { ten_prop_to_write(i, 0, 0) = proprty_col(i); } ten_prop_to_write = ten_prop_to_write.reshape( - Eigen::array{sp_mesh_->GetLx(), - sp_mesh_->GetLy(), - sp_mesh_->GetLz()}); + Eigen::array < Eigen::DenseIndex, 3 > {sp_mesh_->GetLx(), + sp_mesh_->GetLy(), + sp_mesh_->GetLz()}); return ten_prop_to_write; } @@ -339,7 +341,8 @@ namespace da::sha { } } } - std::vector> v_tri = Vec2Triplet(iH, jH, sH); + std::vector > v_tri = Vec2Triplet(iH, jH, + sH); H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles()); H_.setFromTriplets(v_tri.begin(), v_tri.end()); diff --git a/src/Top3d.h b/src/Top3d.h index bfafe95..82a5db3 100644 --- a/src/Top3d.h +++ b/src/Top3d.h @@ -144,7 +144,7 @@ namespace da::sha { } } - Tensor3d TopOptMainLoop(bool only_simulation=false); + Tensor3d TopOptMainLoop(bool only_simulation=false,bool save_internal_rho=false); std::vector @@ -214,6 +214,7 @@ namespace da::sha { public: std::vector v_compliance_; std::vector v_volume_; + std::vector v_ten_rho_; private: std::shared_ptr sp_para_; std::shared_ptr sp_fea_;