diff --git a/README.md b/README.md
index cee79a2..a9f9601 100644
--- a/README.md
+++ b/README.md
@@ -70,6 +70,17 @@ 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 sequence vtk in Paraview.
+4.
### 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 <sp_mech_top3d->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 <sstream>
#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<ConfigMechanicalInterface> sp_mech_inter = std::make_shared<ConfigMechanicalInterface>(
- config_file,tr_model);
- std::shared_ptr<ConfigThermalInterface> sp_ther_inter = std::make_shared<ConfigThermalInterface>(
- config_file,tr_model);
+ std::shared_ptr <ConfigMechanicalInterface> sp_mech_inter = std::make_shared<ConfigMechanicalInterface>(
+ config_file, tr_model);
+ std::shared_ptr <ConfigThermalInterface> sp_ther_inter = std::make_shared<ConfigThermalInterface>(
+ 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<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;
@@ -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<Eigen::DenseIndex, 3>{
- 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<Eigen::DenseIndex, 3>{
- 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<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;
@@ -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<Eigen::DenseIndex, 3>{
- 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<Eigen::DenseIndex, 3>{
- 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<double> v_compliance_;
std::vector<double> v_volume_;
+ std::vector<Tensor3d> v_ten_rho_;
private:
std::shared_ptr<Top3d> 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<Eigen::DenseIndex, 3>{
- 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<Eigen::DenseIndex, 3>{
- 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<Tensor3d>
+ std::vector <Tensor3d>
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<Tensor3d> vt;
+ std::vector <Tensor3d> 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<Eigen::MatrixXd> v_B(8);
+ std::vector <Eigen::MatrixXd> 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;i<delta_coord.rows();++i){
- v_B[i]=sp_fea_->computeBe(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;j<node_id_in_ele_i.size();++j){
- int i_node=node_id_in_ele_i(j);
- if(node_to_write(i_node)!=0){
+ Eigen::VectorXi node_id_in_ele_i = sp_mesh_->MapEleId2NodeIds(
+ 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<Eigen::DenseIndex, 3>{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<Eigen::Triplet<double>> v_tri = Vec2Triplet(iH, jH, sH);
+ std::vector <Eigen::Triplet<double>> 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<Eigen::DenseIndex, 3>{
- 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<Eigen::DenseIndex, 3>{
- 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<Tensor3d>
+ std::vector <Tensor3d>
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<Tensor3d> vt;
+ std::vector <Tensor3d> 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<Eigen::MatrixXd> v_B(8);
+ std::vector <Eigen::MatrixXd> 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;i<delta_coord.rows();++i){
- v_B[i]=sp_fea_->computeBe(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;j<node_id_in_ele_i.size();++j){
- int i_node=node_id_in_ele_i(j);
- if(node_to_write(i_node)!=0){
+ Eigen::VectorXi node_id_in_ele_i = sp_mesh_->MapEleId2NodeIds(
+ 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<Eigen::DenseIndex, 3>{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<Eigen::Triplet<double>> v_tri = Vec2Triplet(iH, jH, sH);
+ std::vector <Eigen::Triplet<double>> 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<Tensor3d>
@@ -214,6 +214,7 @@ namespace da::sha {
public:
std::vector<double> v_compliance_;
std::vector<double> v_volume_;
+ std::vector<Tensor3d> v_ten_rho_;
private:
std::shared_ptr<CtrlPara> sp_para_;
std::shared_ptr<FEA> sp_fea_;