cflin
/
Thermo-elasticTopologyOptimization
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

before debug2

multiple_top
cflin 2 years ago
parent
cafe7ba08d
commit
f5a826a54c
3 changed files with 33 additions and 55 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 16
      examples/top-thermoelastic-regular-model/config.json
  2. 67
      src/ThermoelasticTop3d.cpp
  3. 5
      src/Top3d.cpp

16
examples/top-thermoelastic-regular-model/config.json
View File

@ -11,10 +11,10 @@
"penal": 3.0, "penal": 3.0,
"r_min": 1.5, "r_min": 1.5,
"T_ref": 295, "T_ref": 295,
"T_limit": 345, "T_limit": 325,
"R_E": 12, "R_E": 8,
"R_lambda": 8, "R_lambda": 8,
"R_beta": 0 "R_beta":0
}, },
"model": { "model": {
"regular_model": { "regular_model": {
@ -29,7 +29,7 @@
{ {
"min": [0, 0.5, 0.5], "min": [0, 0.5, 0.5],
"max": [1, 0.5, 0.5], "max": [1, 0.5, 0.5],
"val": [0.0, 0.0, -1e15] "val": [0.0, 0.0, -1.5e9]
} }
], ],
@ -51,12 +51,12 @@
{ {
"min": [0, 0.2, 0.5], "min": [0, 0.2, 0.5],
"max": [1, 0.2, 0.5], "max": [1, 0.2, 0.5],
"heat_flux": 0.335, "unit": "W" "heat_flux": 12.35, "unit": "W"
}, },
{ {
"min": [0, 0.8, 0.5], "min": [0, 0.8, 0.5],
"max": [1, 0.8, 0.5], "max": [1, 0.8, 0.5],
"heat_flux": 0.335, "unit": "W" "heat_flux": 12.35, "unit": "W"
} }
], ],
@ -64,12 +64,12 @@
{ {
"min": [0, 0, 0], "min": [0, 0, 0],
"max": [1, 1, 0], "max": [1, 1, 0],
"temperature": 325, "unit": "K" "temperature": 312, "unit": "K"
}, },
{ {
"min": [0, 0, 1], "min": [0, 0, 1],
"max": [1, 1, 1], "max": [1, 1, 1],
"temperature": 325, "unit": "K" "temperature": 312, "unit": "K"
} }
] ]
} }

67
src/ThermoelasticTop3d.cpp
View File

@ -34,7 +34,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
Eigen::VectorXd dv(sp_mesh_->GetNumEles()); Eigen::VectorXd dv(sp_mesh_->GetNumEles());
dv.setOnes(); dv.setOnes();
dv = sp_mech_top3d_->H_ * (dv.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); dv = (sp_mech_top3d_->H_ * dv).array() / sp_mech_top3d_->Hs_.array();
Eigen::VectorXd ele_to_write = Eigen::VectorXd ele_to_write =
Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz());
@ -77,10 +77,10 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
spdlog::warn("using Eigen built-in direct solver!"); spdlog::warn("using Eigen built-in direct solver!");
#endif #endif
// start iteration // start iteration
while (change > sp_para_->tol_x * 1e-6 && loop < sp_para_->max_loop) { while (change > sp_para_->tol_x && loop < sp_para_->max_loop) {
++loop; ++loop;
// filter // filter
xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); xPhys_col = (sp_mech_top3d_->H_ * xPhys_col).array() / sp_mech_top3d_->Hs_.array();
auto CalR = [](double rho, double R) { auto CalR = [](double rho, double R) {
return rho / (1.0 + R * (1.0 - rho)); return rho / (1.0 + R * (1.0 - rho));
}; };
@ -93,7 +93,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
}; };
auto CalDRDrho_Vec = [](const Eigen::VectorXd &vec_rho, double R) -> Eigen::VectorXd { auto CalDRDrho_Vec = [](const Eigen::VectorXd &vec_rho, double R) -> Eigen::VectorXd {
auto down = 1 + R * (1 - vec_rho.array()); auto down = 1 + R * (1 - vec_rho.array());
return (1 + R) / down.pow(2); return (1 + R) / (down * down);
}; };
auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd {
return E_min + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); return E_min + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min);
@ -154,7 +154,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
double beta_rho = CalBeta(xPhys_col(i)); double beta_rho = CalBeta(xPhys_col(i));
F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;
} }
if (loop ==1 || loop == sp_para_->max_loop) if (loop == 1 || loop == sp_para_->max_loop)
spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm()); spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm());
Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th;
@ -169,9 +169,9 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
// compliance // compliance
Eigen::VectorXd ce(sp_mesh_->GetNumEles()); Eigen::VectorXd ce(sp_mech_top3d_->sp_mesh_->GetNumEles());
for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { for (int i = 0; i < sp_mech_top3d_->sp_mesh_->GetNumEles(); ++i) {
Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); Eigen::VectorXi dofs_in_ele_i = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXd Ue = sp_mech_top3d_->U_(dofs_in_ele_i); Eigen::VectorXd Ue = sp_mech_top3d_->U_(dofs_in_ele_i);
ce(i) = Ue.transpose() * sp_mech_top3d_->Ke_ * Ue; ce(i) = Ue.transpose() * sp_mech_top3d_->Ke_ * Ue;
} }
@ -190,13 +190,12 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
Eigen::VectorXd v_dFth_drho(i_dFth_drho_.size()); Eigen::VectorXd v_dFth_drho(i_dFth_drho_.size());
for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) {
Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); // Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i);
double Te = T(dofs_th).mean(); double Te = T(dofs_th).mean();
// double beta_rho = CalBeta(xPhys_col(i)); // double beta_rho = CalBeta(xPhys_col(i));
// F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; // F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;
Eigen::VectorXd ele_dFth_drho = Eigen::VectorXd ele_dFth_drho =
(CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;// 24x1
) * Inted_;// 24x1
assert(ele_dFth_drho.size() == 24); assert(ele_dFth_drho.size() == 24);
v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho; v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho;
} }
@ -209,27 +208,27 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
sp_mech_top3d_->sp_mesh_->GetNumDofs()); sp_mech_top3d_->sp_mesh_->GetNumDofs());
Eigen::VectorXd v_dFth_dT(i_dFth_dT_.size()); Eigen::VectorXd v_dFth_dT(i_dFth_dT_.size());
for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) {
Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); // Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); // Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i);
double beta_rho = CalBeta(xPhys_col(i)); double beta_rho = CalBeta(xPhys_col(i));
// F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; // F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;
Eigen::MatrixXd ele_dFth_dT = Eigen::MatrixXd ele_dFth_dT =
Eigen::VectorXd::Ones(dofs_th.size()) * 1.0 / 8.0 * beta_rho * Inted_.transpose(); Eigen::VectorXd::Ones(sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE()) * 1.0 /
sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE() * beta_rho * Inted_.transpose();
assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24);
v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), ele_dFth_dT.size())) = ele_dFth_dT.reshaped(); v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), ele_dFth_dT.size())) = ele_dFth_dT.reshaped();
} }
auto v_dFth_dT_tri = Vec2Triplet(i_dFth_dT_, j_dFth_dT_, v_dFth_dT); auto v_dFth_dT_tri = Vec2Triplet(i_dFth_dT_, j_dFth_dT_, v_dFth_dT);
dFth_dT.setFromTriplets(v_dFth_dT_tri.begin(), v_dFth_dT_tri.end()); dFth_dT.setFromTriplets(v_dFth_dT_tri.begin(), v_dFth_dT_tri.end());
Eigen::VectorXd rhs = dFth_dT *2* lambda_m; Eigen::VectorXd rhs = dFth_dT * 2 * lambda_m;
for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) { for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) {
auto [dof, value] = dof_value; auto [dof, value] = dof_value;
rhs(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value; rhs(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value;
} }
// lambda_t // lambda_t
Eigen::VectorXd lambda_t = solver_th.solve(rhs); Eigen::VectorXd lambda_t = solver_th.solve(rhs);
// dF_drho
Eigen::SparseMatrix<double> &dF_drho = dFth_drho;
// lambda_m_Mul_dKm_drho_Mul_U // lambda_m_Mul_dKm_drho_Mul_U
Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U = Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U =
-CalDEDrho_Vec(xPhys_col).array() * ce.array(); -CalDEDrho_Vec(xPhys_col).array() * ce.array();
@ -251,8 +250,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
+ +
lambda_m_Mul_dKm_drho_Mul_U lambda_m_Mul_dKm_drho_Mul_U
+ +
2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_) 2 * Eigen::VectorXd(dFth_drho * sp_mech_top3d_->U_);
;
// dT_drho // dT_drho
auto CalDKth_Mul_T__For_DTi_Drhoj = [&](int rho_i) -> Eigen::SparseVector<double> { auto CalDKth_Mul_T__For_DTi_Drhoj = [&](int rho_i) -> Eigen::SparseVector<double> {
Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(rho_i); Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(rho_i);
@ -266,7 +264,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
return sp_dKth_Mul_T; return sp_dKth_Mul_T;
}; };
auto CalDTi_Drhoj = [&](int Tdof, const Eigen::SparseVector<double> &sp_dKth_Mul_T) -> double { auto CalDTi_Drhoj = [&](int Tdof, const Eigen::SparseVector<double> &sp_dKth_Mul_T) -> double {
Eigen::VectorXd Li=Eigen::VectorXd ::Zero(sp_dKth_Mul_T.rows()); Eigen::VectorXd Li = Eigen::VectorXd::Zero(sp_dKth_Mul_T.rows());
Li(Tdof) = -1; Li(Tdof) = -1;
for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) { for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) {
auto [dof, value] = dof_value; auto [dof, value] = dof_value;
@ -279,38 +277,19 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
sp_thermal_top3d_->set_dofs_to_load.size()); sp_thermal_top3d_->set_dofs_to_load.size());
for (auto it = map_ele2Limit.begin(); it != map_ele2Limit.end(); ++it) { for (auto it = map_ele2Limit.begin(); it != map_ele2Limit.end(); ++it) {
auto [ele_id, v_limited] = *it; auto [ele_id, v_limited] = *it;
auto sp_dKth_Mul_T=CalDKth_Mul_T__For_DTi_Drhoj(ele_id); auto sp_dKth_Mul_T = CalDKth_Mul_T__For_DTi_Drhoj(ele_id);
for (auto &limited: v_limited) { for (auto &limited: v_limited) {
dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj(limited.dof,sp_dKth_Mul_T); dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj(limited.dof, sp_dKth_Mul_T);
} }
} }
// for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) {
// Eigen::VectorXi dofs_in_ele_i = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i);
// Eigen::VectorXd dKe_th_Mul_T =
// CalDlambdaDrho(xPhys_col(i)) * sp_thermal_top3d_->Ke_ * T(dofs_in_ele_i);
// Eigen::MatrixXd Ke_th(sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE(),
// sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE());
// for (int j1 = 0; j1 < Ke_th.cols(); ++j1) {
// for (int i1 = 0; i1 < Ke_th.rows(); ++i1) {
// Ke_th(i1, j1) = sp_thermal_top3d_->K_.coeffRef(i1, j1);
// }
// }
// Eigen::VectorXd ele_dT_drho = Ke_th.llt().solve(-dKe_th_Mul_T);
// dT_drho(i, dofs_in_ele_i) = ele_dT_drho.transpose();
// }
// for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) {
// auto [dof, value] = dof_value;
// dT_drho.col(dof).setZero();
// }
// dc_dx // dc_dx
Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho;
// dT_dx // dT_dx
Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho;
// mma solver // mma solver
#define SENSITIVITY_SCALE_COEF 100 #define SENSITIVITY_SCALE_COEF 200
size_t num_constraints = size_t num_constraints =
1 + dT_dx.cols();// volume and temperature constraints 1 + dT_dx.cols();// volume and temperature constraints
@ -320,8 +299,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col;
double f0val = c; double f0val = c;
Eigen::VectorXd df0dx = flg_chosen Eigen::VectorXd df0dx = flg_chosen
? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff()*SENSITIVITY_SCALE_COEF
: dc_dx / dc_dx.cwiseAbs().maxCoeff() ; : dc_dx / dc_dx.cwiseAbs().maxCoeff()*SENSITIVITY_SCALE_COEF;
// double fval = v - num_variables * sp_para_->volfrac; // double fval = v - num_variables * sp_para_->volfrac;
Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac), Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac),

5
src/Top3d.cpp
View File

@ -187,8 +187,7 @@ namespace da::sha {
void Top3d::Precompute() { void Top3d::Precompute() {
Eigen::MatrixXi mat_ele2dofs = sp_mesh_->GetEleId2DofsMap(); Eigen::MatrixXi mat_ele2dofs = sp_mesh_->GetEleId2DofsMap();
int dofs_each_ele = sp_mesh_->Get_NUM_NODES_EACH_ELE() * int dofs_each_ele = sp_mesh_->Get_DOFS_EACH_ELE();// 24 for mathe; 8 for heat
sp_mesh_->Get_DOFS_EACH_NODE(); // 24 for mathe; 8 for heat
iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele)) iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele))
.transpose() .transpose()
.reshaped(); .reshaped();
@ -249,7 +248,7 @@ namespace da::sha {
Eigen::SparseMatrix<double> sp_inv_Hs(i_Hs.size(), i_Hs.size()); Eigen::SparseMatrix<double> sp_inv_Hs(i_Hs.size(), i_Hs.size());
auto v_inv_Hs_tri = Vec2Triplet(i_Hs, i_Hs, Eigen::VectorXd(1.0 / Hs_.array())); auto v_inv_Hs_tri = Vec2Triplet(i_Hs, i_Hs, Eigen::VectorXd(1.0 / Hs_.array()));
sp_inv_Hs.setFromTriplets(v_inv_Hs_tri.begin(), v_inv_Hs_tri.end()); sp_inv_Hs.setFromTriplets(v_inv_Hs_tri.begin(), v_inv_Hs_tri.end());
drho_dx_ = sp_inv_Hs * H_; drho_dx_ = sp_inv_Hs * H_.transpose();
} }
} // namespace top } // namespace top

Write Preview
Loading…
Cancel
Save