diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index ace787b..2c49339 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -3,3 +3,362 @@ // #include "ThermoelasticTop3d.h" + +da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { + 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(sp_para_->volfrac); + bool flg_chosen = false; + Eigen::VectorXi chosen_ele_id; +// Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); +// bool flg_chosen = chosen_ele_id.size() != 0; +// if (!flg_chosen) { +// // no part chosen +// xPhys_col.setConstant(sp_para_->volfrac); +// } else { +// // pick chosen part to sim +// xPhys_col = sp_mesh_->GetInitEleRho(); +// xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac); +// } + + int loop = 0; + double change = 1.0; + double E0_m = sp_mech_top3d_->sp_fea_->sp_material_->E; + 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; + // Precompute + + Eigen::VectorXd dv(sp_mesh_->GetNumEles()); + dv.setOnes(); + dv = sp_mech_top3d_->H_ * (dv.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); + + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); + Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + + // dofs of limited T + struct LimitedDof { + int dof; + int idx_of_load_dof; + int idx_in_ele; + + LimitedDof(int dof, int idx_of_load_dof, int idx_in_ele) : dof(dof), idx_of_load_dof(idx_of_load_dof), + idx_in_ele(idx_in_ele) {} + }; + std::map> map_ele2Limit; + std::vector v_dof(sp_thermal_top3d_->set_dofs_to_load.begin(), + sp_thermal_top3d_->set_dofs_to_load.end()); + { + Eigen::MatrixXi ele2dof_map = sp_thermal_top3d_->sp_mesh_->GetEleId2DofsMap(); + // loop ele2dof_map + for (int i = 0; i < ele2dof_map.rows(); ++i) { + for (int j = 0; j < ele2dof_map.cols(); ++j) { + for (int k = 0; k < v_dof.size(); ++k) { + if (ele2dof_map(i, j) == v_dof[k]) { + if (map_ele2Limit.find(i) == map_ele2Limit.end()) { + map_ele2Limit[i] = {LimitedDof(v_dof[k], k, j)}; + } else { + map_ele2Limit[i].push_back(LimitedDof(v_dof[k], k, j)); + } + } + } + } + } + } + spdlog::info("end Precompute"); + +#ifdef USE_SUITESPARSE + spdlog::info("using suitesparse solver"); +#else + spdlog::warn("using Eigen built-in direct solver!"); +#endif + // start iteration + while (change > sp_para_->tol_x * 1 && loop < sp_para_->max_loop) { + ++loop; + // filter + xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); + auto CalR = [](double rho, double R) { + return rho / (1.0 + R * (1.0 - rho)); + }; + auto CalR_Vec = [](const Eigen::VectorXd &vec_rho, double R) -> Eigen::VectorXd { + return vec_rho.array() / (1.0 + R * (1.0 - vec_rho.array())); + }; + auto CalDRDrho = [](double rho, double R) { + double down = 1 + R * (1 - rho); + 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()); + return (1 + R) / down.pow(2); + }; + auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + return CalR_Vec(vec_rho, sp_para_->R_E) * E0_m; + }; + auto CalDEDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * E0_m; + }; + auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + return lambda_min + CalR_Vec(vec_rho, sp_para_->R_lambda).array() * (lambda0 - lambda_min); + }; + auto CalDlambdaDrho = [&](double rho) { + 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); + }; + + auto CalBeta = [&](double rho) { + return CalR(rho, sp_para_->R_beta) * E0_m * alpha0; + }; + auto CalDBetaDrho = [&](double rho) { + return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; + }; + // solve T + + Eigen::VectorXd sK_th = + (sp_thermal_top3d_->sKe_ * CalLambda_Vec(xPhys_col).transpose()) + .reshaped(); + auto v_tri_th = Vec2Triplet(sp_thermal_top3d_->iK_, sp_thermal_top3d_->jK_, sK_th); + sp_thermal_top3d_->K_.setFromTriplets(v_tri_th.begin(), v_tri_th.end()); + sp_thermal_top3d_->IntroduceFixedDofs(sp_thermal_top3d_->K_, sp_thermal_top3d_->F_); +#ifdef USE_SUITESPARSE + Eigen::CholmodSupernodalLLT> solver_th; +#else + Eigen::SimplicialLLT> solver_th; +#endif + solver_th.compute(sp_thermal_top3d_->K_); + sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); + + + // solve U + + Eigen::VectorXd sK_m = + (sp_mech_top3d_->sKe_ * CalE_Vec(xPhys_col).transpose()) + .reshaped(); + auto v_tri_m = Vec2Triplet(sp_mech_top3d_->iK_, sp_mech_top3d_->jK_, sK_m); + sp_mech_top3d_->K_.setFromTriplets(v_tri_m.begin(), v_tri_m.end()); + // for each element + Eigen::VectorXd &T = sp_thermal_top3d_->U_; + Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); + + 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); + 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 F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; + sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); +#ifdef USE_SUITESPARSE + Eigen::CholmodSupernodalLLT> solver; +#else + Eigen::SimplicialLLT> solver; +#endif + solver.compute(sp_mech_top3d_->K_); + sp_mech_top3d_->U_ = solver.solve(F); + + + // 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 Ue = sp_mech_top3d_->U_(dofs_in_ele_i); + ce(i) = Ue.transpose() * sp_mech_top3d_->Ke_ * Ue; + } + double c = + ce.transpose() * CalE_Vec(xPhys_col); + double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); + + // sensitivity + // lambda_m + Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; + // dFth_drho + Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), + sp_mech_top3d_->sp_mesh_->GetNumDofs()); + 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); + 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 + assert(ele_dFth_drho.size() == 24); + v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.rows(), ele_dFth_drho.size())) = ele_dFth_drho; + } + auto v_dFth_drho_tri = Vec2Triplet(i_dFth_drho_, j_dFth_drho_, v_dFth_drho); + dFth_drho.setFromTriplets(v_dFth_drho_tri.begin(), v_dFth_drho_tri.end()); + + // dFth_dT + Eigen::SparseMatrix dFth_dT(sp_thermal_top3d_->sp_mesh_->GetNumDofs(), + 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); + 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(); + assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); + v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.rows(), 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 * 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(); + // lambda_t_Mul_dKt_drho_Mul_T + Eigen::VectorXd ce_th(sp_thermal_top3d_->sp_mesh_->GetNumEles()); + 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 Te = sp_thermal_top3d_->U_(dofs_in_ele_i); + Eigen::VectorXd lambda_t_e = lambda_t(dofs_in_ele_i); + ce_th(i) = lambda_t_e.transpose() * sp_thermal_top3d_->Ke_ * Te; + } + Eigen::VectorXd lambda_t_Mul_dKt_drho_Mul_T = CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array(); + // dc_drho + Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + + lambda_m_Mul_dKm_drho_Mul_U + + 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); +// Eigen::VectorXd dc_drho = +// lambda_m_Mul_dKm_drho_Mul_U + +// 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); + // dT_drho + 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); + for (auto &limited: v_limited) { + dT_drho(ele_id, limited.idx_of_load_dof) = ele_dT_drho(limited.idx_in_ele); + } + } +// 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 = drho_dx_ * dc_drho; + // dT_dx + Eigen::MatrixXd dT_dx = drho_dx_ * dT_drho; + + // mma solver + size_t num_constraints = + 1 + dT_dx.cols();// volume and temperature constraints + size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); + auto mma = std::make_shared(num_variables, num_constraints); + 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(); +// 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(); +// 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(); + static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); + static Eigen::VectorXd up_bounds = Eigen::VectorXd::Ones(num_variables); + + // spdlog::info("mma update"); + mma->Update(variables_tmp.data(), df0dx.data(), fval.data(), dfdx.data(), low_bounds.data(), + up_bounds.data()); + if (flg_chosen) { + change = (variables_tmp - xPhys_col(chosen_ele_id)).cwiseAbs().maxCoeff(); + xPhys_col(chosen_ele_id) = variables_tmp; + } else { + change = (variables_tmp - xPhys_col).cwiseAbs().maxCoeff(); + xPhys_col = variables_tmp; + } + + spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); + std::cout << fval.transpose() << std::endl; +#ifdef WRITE_TENSOR_IN_LOOP + // extract vtk + 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{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + top::WriteTensorToVtk( + da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), + ten_xPhys_to_write, sp_mesh_); +#endif + } + // result + sp_mech_top3d_->rho_ = 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{ + 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(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{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + sp_mech_top3d_->rho_field_one_filled_ = ten_xPhys_to_write; + } + + return sp_mech_top3d_->rho_field_zero_filled_; +} diff --git a/src/ThermoelasticTop3d.h b/src/ThermoelasticTop3d.h index 74683e7..7dcf60c 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -15,364 +15,7 @@ namespace da::sha::top { Precompute(); } - Tensor3d TopOptMainLoop() { - 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(sp_para_->volfrac); - bool flg_chosen = false; - Eigen::VectorXi chosen_ele_id; -// Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); -// bool flg_chosen = chosen_ele_id.size() != 0; -// if (!flg_chosen) { -// // no part chosen -// xPhys_col.setConstant(sp_para_->volfrac); -// } else { -// // pick chosen part to sim -// xPhys_col = sp_mesh_->GetInitEleRho(); -// xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac); -// } - - int loop = 0; - double change = 1.0; - double E0_m = sp_mech_top3d_->sp_fea_->sp_material_->E; - 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; - // Precompute - - Eigen::VectorXd dv(sp_mesh_->GetNumEles()); - dv.setOnes(); - dv = sp_mech_top3d_->H_ * (dv.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); - - Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); - Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); - - // dofs of limited T - struct LimitedDof { - int dof; - int idx_of_load_dof; - int idx_in_ele; - - LimitedDof(int dof, int idx_of_load_dof, int idx_in_ele) : dof(dof), idx_of_load_dof(idx_of_load_dof), - idx_in_ele(idx_in_ele) {} - }; - std::map> map_ele2Limit; - std::vector v_dof(sp_thermal_top3d_->set_dofs_to_load.begin(), - sp_thermal_top3d_->set_dofs_to_load.end()); - { - Eigen::MatrixXi ele2dof_map = sp_thermal_top3d_->sp_mesh_->GetEleId2DofsMap(); - // loop ele2dof_map - for (int i = 0; i < ele2dof_map.rows(); ++i) { - for (int j = 0; j < ele2dof_map.cols(); ++j) { - for (int k = 0; k < v_dof.size(); ++k) { - if (ele2dof_map(i, j) == v_dof[k]) { - if (map_ele2Limit.find(i) == map_ele2Limit.end()) { - map_ele2Limit[i] = {LimitedDof(v_dof[k], k, j)}; - } else { - map_ele2Limit[i].push_back(LimitedDof(v_dof[k], k, j)); - } - } - } - } - } - } - spdlog::info("end Precompute"); - -#ifdef USE_SUITESPARSE - spdlog::info("using suitesparse solver"); -#else - spdlog::warn("using Eigen built-in direct solver!"); -#endif - // start iteration - while (change > sp_para_->tol_x * 1 && loop < sp_para_->max_loop) { - ++loop; - // filter - xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); - auto CalR = [](double rho, double R) { - return rho / (1.0 + R * (1.0 - rho)); - }; - auto CalR_Vec = [](const Eigen::VectorXd &vec_rho, double R) -> Eigen::VectorXd { - return vec_rho.array() / (1.0 + R * (1.0 - vec_rho.array())); - }; - auto CalDRDrho = [](double rho, double R) { - double down = 1 + R * (1 - rho); - 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()); - return (1 + R) / down.pow(2); - }; - auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return CalR_Vec(vec_rho, sp_para_->R_E) * E0_m; - }; - auto CalDEDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * E0_m; - }; - auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return lambda_min + CalR_Vec(vec_rho, sp_para_->R_lambda).array() * (lambda0 - lambda_min); - }; - auto CalDlambdaDrho = [&](double rho) { - 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); - }; - - auto CalBeta = [&](double rho) { - return CalR(rho, sp_para_->R_beta) * E0_m * alpha0; - }; - auto CalDBetaDrho = [&](double rho) { - return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; - }; - // solve T - - Eigen::VectorXd sK_th = - (sp_thermal_top3d_->sKe_ * CalLambda_Vec(xPhys_col).transpose()) - .reshaped(); - auto v_tri_th = Vec2Triplet(sp_thermal_top3d_->iK_, sp_thermal_top3d_->jK_, sK_th); - sp_thermal_top3d_->K_.setFromTriplets(v_tri_th.begin(), v_tri_th.end()); - sp_thermal_top3d_->IntroduceFixedDofs(sp_thermal_top3d_->K_, sp_thermal_top3d_->F_); -#ifdef USE_SUITESPARSE - Eigen::CholmodSupernodalLLT> solver_th; -#else - Eigen::SimplicialLLT> solver_th; -#endif - solver_th.compute(sp_thermal_top3d_->K_); - sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); - - - // solve U - - Eigen::VectorXd sK_m = - (sp_mech_top3d_->sKe_ * CalE_Vec(xPhys_col).transpose()) - .reshaped(); - auto v_tri_m = Vec2Triplet(sp_mech_top3d_->iK_, sp_mech_top3d_->jK_, sK_m); - sp_mech_top3d_->K_.setFromTriplets(v_tri_m.begin(), v_tri_m.end()); - // for each element - Eigen::VectorXd &T = sp_thermal_top3d_->U_; - Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); - - 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); - 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 F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; - sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); -#ifdef USE_SUITESPARSE - Eigen::CholmodSupernodalLLT> solver; -#else - Eigen::SimplicialLLT> solver; -#endif - solver.compute(sp_mech_top3d_->K_); - sp_mech_top3d_->U_ = solver.solve(F); - - - // 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 Ue = sp_mech_top3d_->U_(dofs_in_ele_i); - ce(i) = Ue.transpose() * sp_mech_top3d_->Ke_ * Ue; - } - double c = - ce.transpose() * CalE_Vec(xPhys_col); - double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); - - // sensitivity - // lambda_m - Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; - // dFth_drho - Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), - sp_mech_top3d_->sp_mesh_->GetNumDofs()); - 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); - 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 - assert(ele_dFth_drho.size() == 24); - v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.rows(), ele_dFth_drho.size())) = ele_dFth_drho; - } - auto v_dFth_drho_tri = Vec2Triplet(i_dFth_drho_, j_dFth_drho_, v_dFth_drho); - dFth_drho.setFromTriplets(v_dFth_drho_tri.begin(), v_dFth_drho_tri.end()); - - // dFth_dT - Eigen::SparseMatrix dFth_dT(sp_thermal_top3d_->sp_mesh_->GetNumDofs(), - 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); - 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(); - assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); - v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.rows(), 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 * 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(); - // lambda_t_Mul_dKt_drho_Mul_T - Eigen::VectorXd ce_th(sp_thermal_top3d_->sp_mesh_->GetNumEles()); - 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 Te = sp_thermal_top3d_->U_(dofs_in_ele_i); - Eigen::VectorXd lambda_t_e = lambda_t(dofs_in_ele_i); - ce_th(i) = lambda_t_e.transpose() * sp_thermal_top3d_->Ke_ * Te; - } - Eigen::VectorXd lambda_t_Mul_dKt_drho_Mul_T = CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array(); - // dc_drho - Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + - lambda_m_Mul_dKm_drho_Mul_U + - 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); -// Eigen::VectorXd dc_drho = -// lambda_m_Mul_dKm_drho_Mul_U + -// 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); - // dT_drho - 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); - for (auto &limited: v_limited) { - dT_drho(ele_id, limited.idx_of_load_dof) = ele_dT_drho(limited.idx_in_ele); - } - } -// 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 = drho_dx_ * dc_drho; - // dT_dx - Eigen::MatrixXd dT_dx = drho_dx_ * dT_drho; - - // mma solver - size_t num_constraints = - 1 + dT_dx.cols();// volume and temperature constraints - size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); - auto mma = std::make_shared(num_variables, num_constraints); - 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(); -// 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(); -// 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(); - static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); - static Eigen::VectorXd up_bounds = Eigen::VectorXd::Ones(num_variables); - - // spdlog::info("mma update"); - mma->Update(variables_tmp.data(), df0dx.data(), fval.data(), dfdx.data(), low_bounds.data(), - up_bounds.data()); - if (flg_chosen) { - change = (variables_tmp - xPhys_col(chosen_ele_id)).cwiseAbs().maxCoeff(); - xPhys_col(chosen_ele_id) = variables_tmp; - } else { - change = (variables_tmp - xPhys_col).cwiseAbs().maxCoeff(); - xPhys_col = variables_tmp; - } - - spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); - std::cout << fval.transpose() << std::endl; -#ifdef WRITE_TENSOR_IN_LOOP - // extract vtk - 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{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); - top::WriteTensorToVtk( - da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), - ten_xPhys_to_write, sp_mesh_); -#endif - } - // result - sp_mech_top3d_->rho_ = 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{ - 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(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{ - sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); - sp_mech_top3d_->rho_field_one_filled_ = ten_xPhys_to_write; - } - - return sp_mech_top3d_->rho_field_zero_filled_; - } + Tensor3d TopOptMainLoop(); Eigen::VectorXd GetU() const { return sp_thermal_top3d_->GetU();