From 6e84519dc7d94604bbf127d526b1374f7ce83a13 Mon Sep 17 00:00:00 2001 From: cflin Date: Thu, 15 Jun 2023 16:57:04 +0800 Subject: [PATCH] todo: fix Emin --- .../top-mech-regular-model-test/config.json | 2 +- .../config.json | 18 +- src/ThermoelasticTop3d.cpp | 378 +++++++++++++++ src/ThermoelasticTop3d.h | 381 +-------------- src/Top3d.cpp | 451 +++++++++--------- src/Top3d.h | 1 + 6 files changed, 623 insertions(+), 608 deletions(-) diff --git a/examples/top-mech-regular-model-test/config.json b/examples/top-mech-regular-model-test/config.json index 9da75b5..62ab77f 100644 --- a/examples/top-mech-regular-model-test/config.json +++ b/examples/top-mech-regular-model-test/config.json @@ -6,7 +6,7 @@ "topology": { "max_loop": 100, "volfrac": 0.5, - "penal": 3.0, + "penal": 1.0, "r_min": 1.5 }, "model": { diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index e381535..ac1b932 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -2,23 +2,23 @@ "material": { "E": 1.0, "poisson_ratio": 0.3, - "thermal_conductivity":1e11, "unit": "W/K", - "thermal_expansion_coefficient": 1e-15, "unit": "1/K" + "thermal_conductivity":0.004, "unit": "W/K", + "thermal_expansion_coefficient": 1e-6, "unit": "1/K" }, "topology": { "max_loop": 100, "volfrac": 0.5, "r_min": 1.5, "T_ref": 295, - "T_limit": 3250, - "R_E": 8, + "T_limit": 325, + "R_E": 0, "R_lambda": 0, "R_beta": 0 }, "model": { "regular_model": { - "lx": 1, - "ly": 40, + "lx": 2, + "ly": 20, "lz": 30 }, "visual_model":"2-box-refined.obj" @@ -28,7 +28,7 @@ { "min": [0, 1, 0], "max": [1, 1, 0], - "val": [0.0, 0.0, -500.0] + "val": [0.0, 0.0, -1.0] } ], @@ -45,12 +45,12 @@ { "min": [0, 0.25, 0.5], "max": [1, 0.25, 0.5], - "heat_flux": 0.06, "unit": "W" + "heat_flux": 1e-15, "unit": "W" }, { "min": [0, 0.75, 0.5], "max": [1, 0.75, 0.5], - "heat_flux": 0.06, "unit": "W" + "heat_flux": 1e-15, "unit": "W" } ], diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index ace787b..b308ec5 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -3,3 +3,381 @@ // #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 && 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_ * (lambda_min + +// xPhys_col.array() / +// (1.0 + sp_para_->R_lambda * (1.0 - xPhys_col.array())) * +// (lambda0 - lambda_min)).matrix().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_ * + (1e-6+ + xPhys_col.array() / + (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * + E0_m + ).matrix() + .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() * (xPhys_col.array() / + (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * + E0_m).matrix(); + 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 = + -((1 + sp_para_->R_E) / (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())).pow(2)) * 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 = sp_mech_top3d_->drho_dx_ * dc_drho; +// // dT_dx +// Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; + + // mma solver + size_t num_constraints = + 1 ;// 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 fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac)).finished(); + Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v - num_variables *sp_para_->volfrac)).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).finished().transpose(); + Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints) + << dv).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 5bee205..5d95745 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -15,378 +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_ * (lambda_min + - xPhys_col.array() / - (1.0 + sp_para_->R_lambda * (1.0 - xPhys_col.array())) * - (lambda0 - lambda_min)).matrix().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_ * - ( - xPhys_col.array() / - (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * - E0_m).matrix() - .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() * (xPhys_col.array() / - (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * - E0_m).matrix(); - 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 = - -((1 + sp_para_->R_E) / (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())).pow(2)) * 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 = ((1 + (lambda0 - lambda_min)) / - (1.0 + - (lambda0 - lambda_min) * (1.0 - xPhys_col.array())).pow( - 2)) * - 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<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(); @@ -414,13 +43,6 @@ namespace da::sha::top { D0_ = sp_mech_top3d_->sp_fea_->computeD(1.0); const Eigen::MatrixXd Be = sp_mech_top3d_->sp_fea_->computeBe(); Inted_ = Be.transpose() * D0_ * (Eigen::VectorXd(6) << 1, 1, 1, 0, 0, 0).finished(); - Eigen::VectorXi i_Hs = Eigen::VectorXi::LinSpaced(sp_mech_top3d_->Hs_.size(), 0, - sp_mech_top3d_->Hs_.size()); - Eigen::SparseMatrix sp_Hs(i_Hs.size(), i_Hs.size()); - auto v_tri = Vec2Triplet(i_Hs, i_Hs, sp_mech_top3d_->Hs_); - sp_Hs.setFromTriplets(v_tri.begin(), v_tri.end()); - drho_dx_ = sp_Hs * sp_mech_top3d_->H_; - } @@ -431,7 +53,6 @@ namespace da::sha::top { Eigen::VectorXi i_dFth_drho_, j_dFth_drho_; Eigen::VectorXd Inted_; Eigen::MatrixXd D0_; - Eigen::SparseMatrix drho_dx_; }; } diff --git a/src/Top3d.cpp b/src/Top3d.cpp index 581d56b..d4aade0 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -7,235 +7,250 @@ #include #include "Eigen/src/Core/Matrix.h" #include "Util.h" + namespace da::sha { -namespace top { -Tensor3d Top3d::TopOptMainLoop() { - Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); - 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 = sp_fea_->sp_material_->E; - double Emin = E0 * sp_para_->E_factor; - - // precompute - Eigen::VectorXd dv(sp_mesh_->GetNumEles()); - dv.setOnes(); - dv = H_ * (dv.array() / 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(); - spdlog::info("end precompute"); - // // clear output dir - // clear_dir(CMAKE_SOURCE_DIR "/output"); + namespace top { + Tensor3d Top3d::TopOptMainLoop() { + Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); + 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 = sp_fea_->sp_material_->E; + double Emin = E0 * sp_para_->E_factor; + + // precompute + Eigen::VectorXd dv(sp_mesh_->GetNumEles()); + dv.setOnes(); + dv = H_ * (dv.array() / 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(); + spdlog::info("end precompute"); + // // clear output dir + // clear_dir(CMAKE_SOURCE_DIR "/output"); #ifdef USE_SUITESPARSE - spdlog::info("using suitesparse solver"); + spdlog::info("using suitesparse solver"); #else - spdlog::warn("using Eigen built-in direct solver!"); + spdlog::warn("using Eigen built-in direct solver!"); #endif - // start iteration - while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { - ++loop; - Eigen::VectorXd sK = - (sKe_ * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix().transpose()) - .reshaped(); - auto v_tri = Vec2Triplet(iK_, jK_, sK); - K_.setFromTriplets(v_tri.begin(), v_tri.end()); - IntroduceFixedDofs(K_, F_); + // start iteration + while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { + ++loop; + // filter + xPhys_col = H_ * (xPhys_col.array() / Hs_.array()).matrix().eval(); + Eigen::VectorXd sK = + (sKe_ * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix().transpose()) + .reshaped(); + auto v_tri = Vec2Triplet(iK_, jK_, sK); + K_.setFromTriplets(v_tri.begin(), v_tri.end()); + IntroduceFixedDofs(K_, F_); #ifdef USE_SUITESPARSE - Eigen::CholmodSupernodalLLT> solver; + Eigen::CholmodSupernodalLLT> solver; #else - Eigen::SimplicialLLT> solver; + Eigen::SimplicialLLT> solver; #endif - solver.compute(K_); - 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 = U_(dofs_in_ele_i); - ce(i) = Ue.transpose() * Ke_ * Ue; - } - double c = - ce.transpose() * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix(); - double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); - - Eigen::VectorXd dc = - -sp_para_->penal * (E0 - Emin) * xPhys_col.array().pow(sp_para_->penal - 1.0) * ce.array(); - - // mma solver - size_t num_constrants = 1; - size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); - auto mma = std::make_shared(num_variables, num_constrants); - Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; - double f0val = c; - Eigen::VectorXd df0dx = flg_chosen - ? dc(chosen_ele_id).eval() / dc(chosen_ele_id).cwiseAbs().maxCoeff() - : dc / dc.cwiseAbs().maxCoeff(); - double fval = v - num_variables * sp_para_->volfrac; - Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; - 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, 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); + solver.compute(K_); + 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 = U_(dofs_in_ele_i); + ce(i) = Ue.transpose() * Ke_ * Ue; + } + double c = + ce.transpose() * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix(); + double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); + + // sensitivity + Eigen::VectorXd dc_drho = + -sp_para_->penal * (E0 - Emin) * xPhys_col.array().pow(sp_para_->penal - 1.0) * ce.array(); + Eigen::VectorXd dc_dx = drho_dx_ * dc_drho; + + // mma solver + size_t num_constrants = 1; + size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); + auto mma = std::make_shared(num_variables, num_constrants); + 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 dfdx = flg_chosen ? dv(chosen_ele_id) : dv; + 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, 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); #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_); + // 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 - 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()}); - 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()}); - rho_field_one_filled_ = ten_xPhys_to_write; - } - - return rho_field_zero_filled_; -} - -std::vector Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) { - Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); - Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); - // stress - Eigen::MatrixXd mat_stress(sp_mesh_->GetNumEles(), 6); - Eigen::MatrixXd B=sp_fea_->computeBe({0, 0, 0}); - for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); - Eigen::VectorXd Ue = U_(dofs_in_ele_i); - mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue; - } - // fill - std::vector vt; - for (int i = 0; i < which_col_of_stress.size(); ++i) { - ele_to_write(pixel_idx) = mat_stress.col(which_col_of_stress(i)); - vt.push_back(GetTensorFromCol(ele_to_write)); - } - return vt; -} - -Tensor3d Top3d::GetTensorFromCol(const Eigen::VectorXd &proprty_col) { - Tensor3d ten_prop_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1); - assert(proprty_col.size() == ten_prop_to_write.size()); - for (int i = 0; i < proprty_col.size(); ++i) { - ten_prop_to_write(i, 0, 0) = proprty_col(i); - } - ten_prop_to_write = ten_prop_to_write.reshape( - Eigen::array{sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); - return ten_prop_to_write; -} - -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 - iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele)) - .transpose() - .reshaped(); - jK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::RowVectorXi::Ones(dofs_each_ele)) - .transpose() - .reshaped(); - Ke_=sp_fea_->computeKe(1.0); - sKe_ = Ke_.reshaped(); - - // precompute filter - Eigen::VectorXi iH = Eigen::VectorXi::Ones( - sp_mesh_->GetNumEles() * std::pow(2.0 * (std::ceil(sp_para_->r_min) - 1.0) + 1, 3)); - Eigen::VectorXi jH = iH; - Eigen::VectorXd sH(iH.size()); - sH.setZero(); - int cnt = 0; - Hs_ = Eigen::VectorXd(sp_mesh_->GetNumEles()); - - int delta = std::ceil(sp_para_->r_min) - 1; - for (int k = 0; k < sp_mesh_->GetLz(); ++k) { - for (int j = 0; j < sp_mesh_->GetLy(); ++j) { - for (int i = 0; i < sp_mesh_->GetLx(); ++i) { - int ele_id0 = sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i, j, k).finished())(0); - if (ele_id0 == -1) { - continue; + } + // result + 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()}); + 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()}); + rho_field_one_filled_ = ten_xPhys_to_write; + } + + return rho_field_zero_filled_; } - for (int k2 = std::max(k - delta, 0); k2 <= std::min(k + delta, sp_mesh_->GetLz() - 1); - ++k2) { - for (int j2 = std::max(j - delta, 0); j2 <= std::min(j + delta, sp_mesh_->GetLy() - 1); - ++j2) { - for (int i2 = std::max(i - delta, 0); i2 <= std::min(i + delta, sp_mesh_->GetLx() - 1); - ++i2) { - int ele_id1 = - sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i2, j2, k2).finished())(0); - if (ele_id1 == -1) { - continue; - } - iH(cnt) = ele_id0; - jH(cnt) = ele_id1; - sH(cnt) = - std::max(0.0, sp_para_->r_min - Eigen::Vector3d(i - i2, j - j2, k - k2).norm()); - Hs_(ele_id0) += sH(cnt); - ++cnt; + + std::vector Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) { + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); + Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + // stress + Eigen::MatrixXd mat_stress(sp_mesh_->GetNumEles(), 6); + Eigen::MatrixXd B = sp_fea_->computeBe({0, 0, 0}); + for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXd Ue = U_(dofs_in_ele_i); + mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue; } - } + // fill + std::vector vt; + for (int i = 0; i < which_col_of_stress.size(); ++i) { + ele_to_write(pixel_idx) = mat_stress.col(which_col_of_stress(i)); + vt.push_back(GetTensorFromCol(ele_to_write)); + } + return vt; + } + + Tensor3d Top3d::GetTensorFromCol(const Eigen::VectorXd &proprty_col) { + Tensor3d ten_prop_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1); + assert(proprty_col.size() == ten_prop_to_write.size()); + for (int i = 0; i < proprty_col.size(); ++i) { + ten_prop_to_write(i, 0, 0) = proprty_col(i); + } + ten_prop_to_write = ten_prop_to_write.reshape( + Eigen::array{sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + return ten_prop_to_write; + } + + 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 + iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele)) + .transpose() + .reshaped(); + jK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::RowVectorXi::Ones(dofs_each_ele)) + .transpose() + .reshaped(); + Ke_ = sp_fea_->computeKe(1.0); + sKe_ = Ke_.reshaped(); + + // precompute filter + Eigen::VectorXi iH = Eigen::VectorXi::Ones( + sp_mesh_->GetNumEles() * std::pow(2.0 * (std::ceil(sp_para_->r_min) - 1.0) + 1, 3)); + Eigen::VectorXi jH = iH; + Eigen::VectorXd sH(iH.size()); + sH.setZero(); + int cnt = 0; + Hs_ = Eigen::VectorXd(sp_mesh_->GetNumEles()); + + int delta = std::ceil(sp_para_->r_min) - 1; + for (int k = 0; k < sp_mesh_->GetLz(); ++k) { + for (int j = 0; j < sp_mesh_->GetLy(); ++j) { + for (int i = 0; i < sp_mesh_->GetLx(); ++i) { + int ele_id0 = sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i, j, k).finished())(0); + if (ele_id0 == -1) { + continue; + } + for (int k2 = std::max(k - delta, 0); k2 <= std::min(k + delta, sp_mesh_->GetLz() - 1); + ++k2) { + for (int j2 = std::max(j - delta, 0); j2 <= std::min(j + delta, sp_mesh_->GetLy() - 1); + ++j2) { + for (int i2 = std::max(i - delta, 0); i2 <= std::min(i + delta, sp_mesh_->GetLx() - 1); + ++i2) { + int ele_id1 = + sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i2, j2, k2).finished())( + 0); + if (ele_id1 == -1) { + continue; + } + iH(cnt) = ele_id0; + jH(cnt) = ele_id1; + sH(cnt) = + std::max(0.0, + sp_para_->r_min - Eigen::Vector3d(i - i2, j - j2, k - k2).norm()); + Hs_(ele_id0) += sH(cnt); + ++cnt; + } + } + } + } + } + } + std::vector> v_tri = Vec2Triplet(iH, jH, sH); + H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles()); + H_.setFromTriplets(v_tri.begin(), v_tri.end()); + + Eigen::VectorXi i_Hs = Eigen::VectorXi::LinSpaced(Hs_.size(), 0, + Hs_.size()); + 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_; + } - } - } - } - std::vector> v_tri = Vec2Triplet(iH, jH, sH); - H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles()); - H_.setFromTriplets(v_tri.begin(), v_tri.end()); -} -} // namespace top + } // namespace top } // namespace da::sha \ No newline at end of file diff --git a/src/Top3d.h b/src/Top3d.h index 30e5e86..017b411 100644 --- a/src/Top3d.h +++ b/src/Top3d.h @@ -188,6 +188,7 @@ namespace da::sha { SpMat H_; Eigen::VectorXd Hs_; + Eigen::SparseMatrix drho_dx_; // result Eigen::VectorXd U_; Eigen::VectorXd rho_;