From f5a826a54cb5110000aa279427e34de496011c15 Mon Sep 17 00:00:00 2001 From: cflin Date: Sun, 25 Jun 2023 18:53:52 +0800 Subject: [PATCH] before debug2 --- .../config.json | 16 ++--- src/ThermoelasticTop3d.cpp | 67 +++++++------------ src/Top3d.cpp | 5 +- 3 files changed, 33 insertions(+), 55 deletions(-) diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index 2d5fd44..9885fa3 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -11,10 +11,10 @@ "penal": 3.0, "r_min": 1.5, "T_ref": 295, - "T_limit": 345, - "R_E": 12, + "T_limit": 325, + "R_E": 8, "R_lambda": 8, - "R_beta": 0 + "R_beta":0 }, "model": { "regular_model": { @@ -29,7 +29,7 @@ { "min": [0, 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], "max": [1, 0.2, 0.5], - "heat_flux": 0.335, "unit": "W" + "heat_flux": 12.35, "unit": "W" }, { "min": [0, 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], "max": [1, 1, 0], - "temperature": 325, "unit": "K" + "temperature": 312, "unit": "K" }, { "min": [0, 0, 1], "max": [1, 1, 1], - "temperature": 325, "unit": "K" + "temperature": 312, "unit": "K" } ] } diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index c79ac45..dcebab2 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -34,7 +34,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd dv(sp_mesh_->GetNumEles()); 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::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!"); #endif // 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; // 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) { 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 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 { 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)); 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()); 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 - Eigen::VectorXd ce(sp_mesh_->GetNumEles()); - for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXd ce(sp_mech_top3d_->sp_mesh_->GetNumEles()); + for (int i = 0; i < sp_mech_top3d_->sp_mesh_->GetNumEles(); ++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); 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()); 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_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 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; } @@ -209,27 +208,27 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { sp_mech_top3d_->sp_mesh_->GetNumDofs()); Eigen::VectorXd v_dFth_dT(i_dFth_dT_.size()); 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_m = sp_mech_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); double beta_rho = CalBeta(xPhys_col(i)); // F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; 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); 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); 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) { auto [dof, value] = dof_value; rhs(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value; } // lambda_t Eigen::VectorXd lambda_t = solver_th.solve(rhs); - // dF_drho - Eigen::SparseMatrix &dF_drho = dFth_drho; + // lambda_m_Mul_dKm_drho_Mul_U Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U = -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 + - 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_) - ; + 2 * Eigen::VectorXd(dFth_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); @@ -266,7 +264,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 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()); + 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; @@ -279,38 +277,19 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 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; - 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) { - 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 Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; // dT_dx Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; // mma solver -#define SENSITIVITY_SCALE_COEF 100 +#define SENSITIVITY_SCALE_COEF 200 size_t num_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; double f0val = c; Eigen::VectorXd df0dx = flg_chosen - ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() - : dc_dx / dc_dx.cwiseAbs().maxCoeff() ; + ? 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; // double fval = v - num_variables * sp_para_->volfrac; Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac), diff --git a/src/Top3d.cpp b/src/Top3d.cpp index d4aade0..cfd3e90 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -187,8 +187,7 @@ namespace da::sha { void Top3d::Precompute() { Eigen::MatrixXi mat_ele2dofs = sp_mesh_->GetEleId2DofsMap(); - int dofs_each_ele = sp_mesh_->Get_NUM_NODES_EACH_ELE() * - sp_mesh_->Get_DOFS_EACH_NODE(); // 24 for mathe; 8 for heat + int dofs_each_ele = sp_mesh_->Get_DOFS_EACH_ELE();// 24 for mathe; 8 for heat iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele)) .transpose() .reshaped(); @@ -249,7 +248,7 @@ namespace da::sha { Eigen::SparseMatrix 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())); 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