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

fixed bugs

multiple_top
cflin 2 years ago
parent
d2787f90df
commit
cafe7ba08d
4 changed files with 49 additions and 32 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      CMakeLists.txt
  2. 5
      examples/top-thermoelastic-regular-model/config.json
  3. 1
      examples/top-thermoelastic-regular-model/main.cpp
  4. 73
      src/ThermoelasticTop3d.cpp

2
CMakeLists.txt
View File

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.25)
cmake_minimum_required(VERSION 3.24)
project(top3d)
option(PROJECT_WITH_SIMD "Enable SIMD" ON)
option(ENABLE_SUITESPARSE "Use SuiteSparse" ON)

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

@ -6,8 +6,9 @@
"thermal_expansion_coefficient": 1.21e-5, "unit": "1/K"
},
"topology": {
"max_loop": 100,
"max_loop": 150,
"volfrac": 0.5,
"penal": 3.0,
"r_min": 1.5,
"T_ref": 295,
"T_limit": 345,
@ -28,7 +29,7 @@
{
"min": [0, 0.5, 0.5],
"max": [1, 0.5, 0.5],
"val": [0.0, 0.0, -1e1]
"val": [0.0, 0.0, -1e15]
}
],

1
examples/top-thermoelastic-regular-model/main.cpp
View File

@ -34,6 +34,7 @@ int main() {
para->max_loop = j_config["topology"]["max_loop"];
para->volfrac = j_config["topology"]["volfrac"];
para->r_min = j_config["topology"]["r_min"];
para->penal=j_config["topology"]["penal"];
para->T_ref=j_config["topology"]["T_ref"];
para->T_limit=j_config["topology"]["T_limit"];
para->R_E=j_config["topology"]["R_E"];

73
src/ThermoelasticTop3d.cpp
View File

@ -29,7 +29,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
double lambda0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_conductivity;
double lambda_min = lambda0 * sp_mech_top3d_->sp_para_->E_factor;
double alpha0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_expansion_coefficient;
double alpha_min=alpha0 * sp_mech_top3d_->sp_para_->E_factor;
double alpha_min = alpha0 * sp_mech_top3d_->sp_para_->E_factor;
// Precompute
Eigen::VectorXd dv(sp_mesh_->GetNumEles());
@ -77,7 +77,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
spdlog::warn("using Eigen built-in direct solver!");
#endif
// start iteration
while (change > sp_para_->tol_x*1e-6 && loop < sp_para_->max_loop) {
while (change > sp_para_->tol_x * 1e-6 && loop < sp_para_->max_loop) {
++loop;
// filter
xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval();
@ -89,7 +89,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
};
auto CalDRDrho = [](double rho, double R) {
double down = 1 + R * (1 - rho);
return (1 + R) / down * down;
return (1 + R) / (down * down);
};
auto CalDRDrho_Vec = [](const Eigen::VectorXd &vec_rho, double R) -> Eigen::VectorXd {
auto down = 1 + R * (1 - vec_rho.array());
@ -108,15 +108,17 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
return CalDRDrho(rho, sp_para_->R_lambda) * (lambda0 - lambda_min);
};
auto CalDlambdaDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd {
return CalDRDrho_Vec(vec_rho, sp_para_->R_lambda) * (lambda0 - lambda_min);
return CalDRDrho_Vec(vec_rho, sp_para_->R_lambda) * (lambda0 - lambda_min);
};
auto CalBeta = [&](double rho) {
return E_min*alpha_min+ CalR(rho, sp_para_->R_beta) * (E0_m * alpha0-E_min*alpha_min);
return CalR(rho, sp_para_->R_beta) * E0_m * alpha0;
};
auto CalDBetaDrho = [&](double rho) {
return CalDRDrho(rho, sp_para_->R_beta) * (E0_m * alpha0-E_min*alpha_min);
return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0;
};
// solve T
Eigen::VectorXd sK_th =
@ -152,7 +154,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
double beta_rho = CalBeta(xPhys_col(i));
F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;
}
spdlog::info("||Fth|| / ||Fm||: {}",F_th.norm()/sp_mech_top3d_->F_.norm());
if (loop ==1 || loop == sp_para_->max_loop)
spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm());
Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th;
sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F);
@ -189,10 +192,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i);
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_;
Eigen::VectorXd ele_dFth_drho =
CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;// 24x1
(CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref)
) * Inted_;// 24x1
assert(ele_dFth_drho.size() == 24);
v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho;
}
@ -217,7 +221,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
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());
Eigen::VectorXd rhs = dFth_dT * lambda_m;
Eigen::VectorXd rhs = dFth_dT *2* lambda_m;
for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) {
auto [dof, value] = dof_value;
rhs(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value;
@ -246,28 +250,38 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
lambda_t_Mul_dKt_drho_Mul_T
+
lambda_m_Mul_dKm_drho_Mul_U
+
2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_)
;
+
2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_)
;
// dT_drho
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::VectorXd ele_T = sp_thermal_top3d_->U_(dofs_in_ele);
Eigen::SparseVector<double> sp_dKth_Mul_T(sp_thermal_top3d_->sp_mesh_->GetNumDofs());
Eigen::VectorXd dKe_th_Mul_T =
CalDlambdaDrho(xPhys_col(rho_i)) * sp_thermal_top3d_->Ke_ * ele_T;
for (int i = 0; i < dofs_in_ele.size(); ++i) {
sp_dKth_Mul_T.coeffRef(dofs_in_ele(i)) = dKe_th_Mul_T(i);
}
return sp_dKth_Mul_T;
};
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());
Li(Tdof) = -1;
for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) {
auto [dof, value] = dof_value;
Li(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value;
}
Eigen::VectorXd lambda_i = solver_th.solve(Li);
return lambda_i.transpose() * sp_dKth_Mul_T;
};
Eigen::MatrixXd dT_drho = Eigen::MatrixXd::Zero(sp_thermal_top3d_->sp_mesh_->GetNumEles(),
sp_thermal_top3d_->set_dofs_to_load.size());
for (auto it = map_ele2Limit.begin(); it != map_ele2Limit.end(); ++it) {
auto [ele_id, v_limited] = *it;
Eigen::VectorXi dofs_in_ele_i = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(ele_id);
Eigen::VectorXd dKe_th_Mul_T =
CalDlambdaDrho(xPhys_col(ele_id)) * 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);
auto sp_dKth_Mul_T=CalDKth_Mul_T__For_DTi_Drhoj(ele_id);
for (auto &limited: v_limited) {
dT_drho(ele_id, limited.idx_of_load_dof) = ele_dT_drho(limited.idx_in_ele);
dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj(limited.dof,sp_dKth_Mul_T);
}
}
@ -306,15 +320,16 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() {
Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col;
double f0val = c;
Eigen::VectorXd df0dx = flg_chosen
? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() *SENSITIVITY_SCALE_COEF
: dc_dx / dc_dx.cwiseAbs().maxCoeff() *SENSITIVITY_SCALE_COEF;
? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff()
: dc_dx / dc_dx.cwiseAbs().maxCoeff() ;
// double fval = v - num_variables * sp_para_->volfrac;
Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac),
T(v_dof).array() / sp_para_->T_limit - 1).finished() * SENSITIVITY_SCALE_COEF;
// Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv;
Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints)
<< 1.0 / num_variables * dv , 1.0 / sp_para_->T_limit * dT_dx).finished().transpose() *SENSITIVITY_SCALE_COEF;
<< 1.0 / num_variables * dv, 1.0 / sp_para_->T_limit * dT_dx).finished().transpose() *
SENSITIVITY_SCALE_COEF;
static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables);

Write Preview
Loading…
Cancel
Save