From 0d5b505c75b4b04ec02b6d9798c2a1784fa4dd08 Mon Sep 17 00:00:00 2001 From: cflin Date: Sat, 1 Jul 2023 19:53:07 +0800 Subject: [PATCH] before add inregular --- CMakeLists.txt | 1 + .../CMakeLists.txt | 2 - .../config.json | 43 ++++------- .../config_bak.json | 76 +++++++++++++++++++ .../top-thermoelastic-regular-model/main.cpp | 8 +- src/ThermoelasticTop3d.cpp | 59 ++++++++++---- src/Top3d.h | 2 +- 7 files changed, 143 insertions(+), 48 deletions(-) create mode 100644 examples/top-thermoelastic-regular-model/config_bak.json diff --git a/CMakeLists.txt b/CMakeLists.txt index 639cc7c..f884178 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -57,5 +57,6 @@ target_link_libraries(${PROJECT_NAME}_lib PUBLIC mma::mma mma::gcmma) target_compile_definitions(${PROJECT_NAME}_lib PUBLIC CMAKE_SOURCE_DIR="${CMAKE_SOURCE_DIR}") target_compile_definitions(${PROJECT_NAME}_lib PUBLIC WRITE_TENSOR) target_compile_definitions(${PROJECT_NAME}_lib PUBLIC DEBUG) +#target_compile_definitions(${PROJECT_NAME}_lib PUBLIC MECH_ONLY) add_subdirectory(${CMAKE_SOURCE_DIR}/examples) \ No newline at end of file diff --git a/examples/top-thermoelastic-regular-model/CMakeLists.txt b/examples/top-thermoelastic-regular-model/CMakeLists.txt index 14567b2..9158481 100644 --- a/examples/top-thermoelastic-regular-model/CMakeLists.txt +++ b/examples/top-thermoelastic-regular-model/CMakeLists.txt @@ -7,5 +7,3 @@ target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC CONFIG_FILE="${CMAKE_CURRE target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC OUTPUT_DIR="${CMAKE_SOURCE_DIR}/output") target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC ASSETS_DIR="${CMAKE_SOURCE_DIR}/assets") target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC DA_CMD) -#target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC ASSETS_DIR="${CMAKE_SOURCE_DIR}/assets") - diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index 9885fa3..c8f1fa6 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -1,19 +1,19 @@ { "material": { - "E": 2e11, + "E": 2.1e11, "poisson_ratio": 0.3, "thermal_conductivity":43, "unit": "W/(m*K)", "thermal_expansion_coefficient": 1.21e-5, "unit": "1/K" }, "topology": { - "max_loop": 150, - "volfrac": 0.5, + "max_loop": 200, + "volfrac": 0.4, "penal": 3.0, - "r_min": 1.5, - "T_ref": 295, - "T_limit": 325, - "R_E": 8, - "R_lambda": 8, + "r_min": 1.1, + "T_ref": 0, + "T_limit": 26, + "R_E": 28, + "R_lambda": 28, "R_beta":0 }, "model": { @@ -27,9 +27,9 @@ "mechanical_boundary_condition":{ "NBC": [ { - "min": [0, 0.5, 0.5], - "max": [1, 0.5, 0.5], - "val": [0.0, 0.0, -1.5e9] + "min": [0, 0.44, 0], + "max": [1, 0.56, 0], + "val": [0.0, 0.0, -8.4e8] } ], @@ -49,27 +49,16 @@ "thermal_boundary_condition": { "NBC": [ { - "min": [0, 0.2, 0.5], - "max": [1, 0.2, 0.5], - "heat_flux": 12.35, "unit": "W" - }, - { - "min": [0, 0.8, 0.5], - "max": [1, 0.8, 0.5], - "heat_flux": 12.35, "unit": "W" + "min": [0, 1, 0.3], + "max": [1, 1, 0.7], + "heat_flux": 5, "unit": "W" } - ], "DBC": [ { "min": [0, 0, 0], - "max": [1, 1, 0], - "temperature": 312, "unit": "K" - }, - { - "min": [0, 0, 1], - "max": [1, 1, 1], - "temperature": 312, "unit": "K" + "max": [1, 0, 1], + "temperature": 0, "unit": "K" } ] } diff --git a/examples/top-thermoelastic-regular-model/config_bak.json b/examples/top-thermoelastic-regular-model/config_bak.json new file mode 100644 index 0000000..c5b6ca8 --- /dev/null +++ b/examples/top-thermoelastic-regular-model/config_bak.json @@ -0,0 +1,76 @@ +{ + "material": { + "E": 2.1e11, + "poisson_ratio": 0.3, + "thermal_conductivity":43, "unit": "W/(m*K)", + "thermal_expansion_coefficient": 1.21e-5, "unit": "1/K" + }, + "topology": { + "max_loop": 200, + "volfrac": 0.5, + "penal": 3.0, + "r_min": 1.05, + "T_ref": 312, + "T_limit": 316.52, + "R_E": 28, + "R_lambda": 28, + "R_beta":0 + }, + "model": { + "regular_model": { + "lx": 1, + "ly": 50, + "lz": 35 + }, + "visual_model":"2-box-refined.obj" + }, + "mechanical_boundary_condition":{ + "NBC": [ + { + "min": [0, 0.5, 0.5], + "max": [1, 0.5, 0.5], + "val": [0.0, 0.0, -4.0e8] + } + + ], + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 0, 1], + "dir": [1, 1, 1] + }, + { + "min": [0, 1, 0], + "max": [1, 1, 1], + "dir": [1, 1, 1] + } + ] + }, + "thermal_boundary_condition": { + "NBC": [ + { + "min": [0, 0.2, 0.5], + "max": [1, 0.2, 0.5], + "heat_flux": 3.35,"back":12.35, "unit": "W" + }, + { + "min": [0, 0.8, 0.5], + "max": [1, 0.8, 0.5], + "heat_flux": 3.35,"back":12.35, "unit": "W" + } + + ], + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 1, 0], + "temperature": 312, "unit": "K" + }, + { + "min": [0, 0, 1], + "max": [1, 1, 1], + "temperature": 312, "unit": "K" + } + ] + } +} \ No newline at end of file diff --git a/examples/top-thermoelastic-regular-model/main.cpp b/examples/top-thermoelastic-regular-model/main.cpp index 8d62a3a..d40e2a2 100644 --- a/examples/top-thermoelastic-regular-model/main.cpp +++ b/examples/top-thermoelastic-regular-model/main.cpp @@ -91,7 +91,8 @@ int main() { Eigen::Vector3d val(NBCI["val"][0], NBCI["val"][1], NBCI["val"][2]); top::Neu t_neu(minBBox, maxBBox, val); - sp_mech_top3d->AddNBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box), t_neu.val); + Eigen::MatrixXi coords=mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_mech_top3d->AddNBC(coords, t_neu.val/coords.rows()); } } @@ -118,8 +119,9 @@ int main() { Eigen::Vector3d maxBBox(NBCI["max"][0], NBCI["max"][1], NBCI["max"][2]); Eigen::Vector3d val( NBCI["heat_flux"],0,0); top::Neu t_neu(minBBox, maxBBox, val); -// auto ttt=mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); - sp_thermal_top3d->AddNBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box), t_neu.val); + + Eigen::MatrixXi coords=mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_thermal_top3d->AddNBC(coords, t_neu.val/coords.rows()); } } diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index dcebab2..bf0eea7 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -8,7 +8,7 @@ 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); + xPhys_col.setConstant(sp_para_->volfrac/sp_para_->volfrac); bool flg_chosen = false; Eigen::VectorXi chosen_ele_id; // Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); @@ -118,7 +118,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; }; - +#ifndef MECH_ONLY // solve T Eigen::VectorXd sK_th = @@ -134,7 +134,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { #endif solver_th.compute(sp_thermal_top3d_->K_); sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); - +#endif // solve U @@ -143,7 +143,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { .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 +#ifndef MECH_ONLY + // F_th Eigen::VectorXd &T = sp_thermal_top3d_->U_; Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); @@ -156,7 +157,13 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { } 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; +#endif + Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + #ifndef MECH_ONLY + + F_th + +#endif + ; sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); #ifdef USE_SUITESPARSE @@ -183,7 +190,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // sensitivity // lambda_m Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; - +#ifndef MECH_ONLY // dFth_drho Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), sp_mech_top3d_->sp_mesh_->GetNumDofs()); @@ -228,10 +235,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { } // lambda_t Eigen::VectorXd lambda_t = solver_th.solve(rhs); - +#endif // lambda_m_Mul_dKm_drho_Mul_U Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U = -CalDEDrho_Vec(xPhys_col).array() * ce.array(); +#ifndef MECH_ONLY // 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) { @@ -245,12 +253,20 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // 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_); +#endif +#ifdef MECH_ONLY + Eigen::VectorXd dc_drho = + lambda_m_Mul_dKm_drho_Mul_U + ; +#else Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + lambda_m_Mul_dKm_drho_Mul_U + 2 * Eigen::VectorXd(dFth_drho * sp_mech_top3d_->U_); +#endif +#ifndef MECH_ONLY // 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); @@ -283,15 +299,20 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { } } - // 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; +#endif + // dc_dx + Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; // mma solver #define SENSITIVITY_SCALE_COEF 200 size_t num_constraints = - 1 + dT_dx.cols();// volume and temperature constraints + 1 + #ifndef MECH_ONLY + + dT_dx.cols()// volume and temperature constraints +#endif + ; size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); @@ -299,15 +320,23 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; double f0val = c; Eigen::VectorXd df0dx = flg_chosen - ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff()*SENSITIVITY_SCALE_COEF - : dc_dx / dc_dx.cwiseAbs().maxCoeff()*SENSITIVITY_SCALE_COEF; + ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() * + 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), - T(v_dof).array() / sp_para_->T_limit - 1).finished() * SENSITIVITY_SCALE_COEF; + Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac) +#ifndef MECH_ONLY + , T(v_dof).array() / sp_para_->T_limit - 1 +#endif + ).finished() * SENSITIVITY_SCALE_COEF; // Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints) - << 1.0 / num_variables * dv, 1.0 / sp_para_->T_limit * dT_dx).finished().transpose() * + << 1.0 / num_variables * dv +#ifndef MECH_ONLY + , 1.0 / sp_para_->T_limit * dT_dx +#endif + ).finished().transpose() * SENSITIVITY_SCALE_COEF; diff --git a/src/Top3d.h b/src/Top3d.h index 017b411..9d9de6f 100644 --- a/src/Top3d.h +++ b/src/Top3d.h @@ -29,7 +29,7 @@ namespace da::sha { double r_min = 2.0; double T_ref = 295; double T_limit=325; - double tol_x = 0.01; + double tol_x = 0.00001; double E_factor = 1e-9; double R_E = 8; double R_lambda = 0;