From cafe7ba08d44ca79c99bf6a449cec6afc5e14460 Mon Sep 17 00:00:00 2001 From: cflin Date: Fri, 16 Jun 2023 03:31:55 +0800 Subject: [PATCH] fixed bugs --- CMakeLists.txt | 2 +- .../config.json | 5 +- .../top-thermoelastic-regular-model/main.cpp | 1 + src/ThermoelasticTop3d.cpp | 73 +++++++++++-------- 4 files changed, 49 insertions(+), 32 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c1101b2..639cc7c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index 5434b63..2d5fd44 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -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] } ], diff --git a/examples/top-thermoelastic-regular-model/main.cpp b/examples/top-thermoelastic-regular-model/main.cpp index 217d795..8d62a3a 100644 --- a/examples/top-thermoelastic-regular-model/main.cpp +++ b/examples/top-thermoelastic-regular-model/main.cpp @@ -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"]; diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index d4af1c9..c79ac45 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -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 { + 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 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 &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);